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INTRODUCTION 


Accuracy  in  the  prediction  of  osteoporotic  fracture  risk  using  bone  mineral  density  alone  is  limited.  Experiments 
using  microcomputed  tomography  (pCT),  a  laboratory  scale  imaging  device,  on  cadaver  bones  show  that 
parameters  that  quantify  the  heterogeneity  of  the  cancellous  bone  microstructure  are  especially  useful  for 
predicting  vertebral  strength.  However,  measurement  of  bone  microstructure  in  a  clinical  setting  is  difficult  and, 
when  possible,  is  generally  limited  to  easily  accessible  extremities.  Digital  tomosynthesis  (DTS)  is  a  clinically 
available  imaging  technique  that  can  produce  multiple  planar  images  of  an  object  with  high  resolution,  as  if  the 
object  was  physically  sliced  before  imaging.  The  objective  of  the  current  project  is  to  determine  the  potential  of 
DTS  for  quantitative  analysis  of  vertebral  bone  quality  and  for  predicting  vertebral  bone  fragility.  As  a  part  of 
this  evaluation,  parameters  derived  from  DTS  are  compared  to  those  measured  from  pCT  and  their  ability  to 
predict  in  vitro  vertebral  strength  is  examined  by  a  set  of  destructive  tests  on  vertebrae  from  human  cadavers. 
Methods  of  fractal,  lacunarity,  mean  intercept  length  and  line  fraction  deviation  analyses  are  used  to  quantify 
cancellous  microstructure  from  DTS  images.  The  strength  of  the  vertebrae  will  be  measured  using  a  loading 
mode  that  will  generate  a  wedge-shaped  fracture  as  seen  in  clinical  vertebral  fractures.  If  DTS  provides 
structural  information  sufficiently  detailed  compared  to  the  pCT  gold  standard  or  that  significantly  correlates 
with  the  strength,  the  development  of  this  concept  into  a  clinical  screening  process  will  be  pursued  through 
further  refining  and  testing. 

BODY 

Task  1.  Retrieval  of  20  cadaveric  human  fresh-frozen  thoraco-lumbar  vertebrae  from  tissue  banks  and  body 
donation  programs  (Months  1-4): 

la.  Review  and  activation  of  tissue  collection  protocols  by  the  donation  program  (Months  1-2). 

lb.  Collection  and  shipment  of  the  vertebrae  (Months  2-4) 

Approval  of  the  institutional  review  board  (IRB  6688)  for  the  project  was  obtained  and  submitted  to  the 
USAMRMC  Office  of  Research  Protections  (ORP),  Human  Research  Protection  Office  (HRPO)  for  review  along 
with  a  claim  of  exemption  from  review  to  use  preexisting  cadaveric  vertebrae.  The  statement  that  the  project  may 
proceed  with  no  further  requirement  for  review  by  the  HRPO  was  received  prior  to  the  beginning  of  the  supported 
period.  In  accordance  with  this  notification  (HRPO  Log  Number  A-16751),  newly  retrieved  vertebrae  as  well  as 
previously  existing  deidentified  cadaveric  vertebrae  are  used. 

We  identified  40  fresh-frozen  thoraco-lumbar  vertebrae  from  10  subjects  previously  collected  through  approved 
protocols  and  allocated  them  to  the  project.  In  addition,  we  established  tissue  collection  protocols  with  the 
National  Disease  Research  Interchange  (NDRI),  the  Anatomical  Donations  Program  of  the  University  of  Michigan 
and  Platinum  Training.  Seven  spines  were  supplied  by  NDRI  and  nine  by  Platinum  Training  (total  16  new  spines). 

Vertebrae  were  separated  from  spines,  stripped  of  surrounding  soft  tissue,  and  posterior  elements  were  resected 
flush  with  the  posterior  aspect  of  the  vertebral  body.  Following  preparation,  vertebral  bodies  were  wrapped  in 
saline-soaked  towels  and  preserved  at  -20  degrees  Celsius  until  scanning  and  mechanical  testing  were 
performed. 

This  task  is  complete. 

Task  2.  Refinement  of  Digital  Tomosynthesis  (DTS)  scan  parameters  (Months  1-4): 

2a.  In  order  to  have  reference  materials  that  resemble  cancellous  bone  but  with  estimated  properties,  prepare 
wire-mesh  references  with  a  wire  thickness  of  50pm  and  varying  mesh  density  with  wire  separation  in  the  order 
of  1mm  and  varying  out-of-plane  separation.  (Month  1) 

2b.  Perform  a  series  of  DTS  scans  on  these  references  using  variations  of  tube  angle,  number  of  projection 
images,  number  of  reconstructed  planes  and  the  incident  exposure  for  each  projection.  (Months  1  -3) 

2c.  Perform  pCT  on  reference  materials.  (Months  1-3) 

2d.  Determine  scan  parameters  that  produce  the  best  agreement  between  pCT  and  DTS.  (Month  4) 

In  order  to  refine  DTS  scanning  parameters,  custom  phantoms  were  created  to  assess  the  accuracy  of  DTS  in 
resolution  of  separation,  thickness,  distances,  in-plane  and  out-of-plane  geometric  distortion,  and  density 
linearity. 
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To  assess  the  minimum  spacing  which  may  be  resolved  by  DTS,  six  phantoms  were  created  comprised  of 
equally  spaced  stacked  squares  of  radiopaque  aluminum  and  radiolucent  acrylic  at  350,  300,  250,  200,  150, 
and  100  pm  nominal  spacing,  as  measured  by  micrometer  during  phantom  preparation  (Figure  la).  The 
squares  were  clamped  and  bound  at  the  edges  using  epoxy  to  create  a  solid  construct.  The  resulting  cross- 
sectional  banding  pattern  thickness  and  frequency  is  comparable  to  human  trabecular  bone.  To  assess  the 
minimum  thickness  which  may  be  resolved  by  DTS,  another  set  of  six  phantoms  was  created  with  thin 
radiopaque  aluminum  layers  (nominally  25,  60,  75,  90,  115  and  120  pm  as  manually  measured  from  pCT 
images)  bound  between  thick  (1.5  mm)  layers  of  radiolucent  acrylic  (Figure  1b). 

Density  linearity  phantoms  were  created  using  serial  dilutions  of  28.2%  (282  mg/ml)  organically  bound  iodine  in 
the  form  of  iothalamate  meglumine,  an  ionic  contrast  agent.  Nominal  concentrations  of  282,  141,  70.5,  35.25, 
17.63,  and  8.81  mg/ml  organically  bound  iodine,  as  well  as  air  and  water,  were  placed  in  a  sealed  well  plate 
(Figure  1c). 

For  demonstration  of  anisotropic  and  edge-enhanced  characteristics  of  DTS,  a  nylon  sphere  (Figure  Id)  and 
cylinder  (Figure  1e)  were  included  in  the  phantom  tray.  To  aid  in  measuring  true  voxel  size  in  the  reconstructed 
volume,  a  geometric  accuracy  phantom  was  created  using  four  1  mm  beads,  placed  in  four  corners  at  35  mm 
apart  (Figure  If).  An  embedded  human  vertebra  was  also  included  in  the  phantom  as  a  realistic  reference 
material  (Figure  1g). 


Figure  1 :  Tray  of  phantoms  to  assess  DTS  resolution,  geometric  distortion,  and  density  linearity. 


An  aluminum  grid  panel  with  hollow  evenly  spaced  cells  with  known  dimensions  (12.7  x  12.7  x  12.7  mm,  Figure 
2)  was  used  to  examine  geometric  distortion  and  uneven  signal  distribution  within  the  reconstructed  DTS 
images. 
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Figure  2:  Aluminum  grid  panel  phantom  to  assess  geometric  distortion  in  x,  y  and  z  directions. 

Density  linearity  phantoms 

In  each  of  5  scanning  sessions,  DTS  scans  and  planar  digital  radiographs  were  performed  on  the  sealed  wells 
containing  serial  dilutions  of  iodine  solution.  As  DTS  scans  were  performed  with  the  phantom  tray  oriented  in  3 
directions,  the  density  phantoms  were  scanned  a  total  of  3  times  per  session,  allowing  assessment  of  session- 
to-session  variability  as  well  as  within-session  variability.  A  single  radiograph  was  taken  in  each  session. 

Images  were  inverted  for  white  on  black  presentation,  and  a  cylindrical  selection  was  drawn  at  the  center  of 
each  well.  Average  gray  value  was  measured  for  each  well  at  a  single  orientation  for  each  session.  Only  the 
central  8  mm  of  each  well  was  analyzed  to  avoid  a  darker  region  of  edge  enhancement  at  the  perimeter  of 
each  1 0  mm  diameter  well,  which  is  due  to  a  high-pass  filter  applied  in  the  direction  of  tomographic  acquisition. 
To  assess  within-session  variability,  measurements  were  repeated  at  three  orientations  per  session. 
Radiographs  were  similarly  analyzed  by  recording  average  gray  values  within  an  8  mm  circle  for  each  well. 

Spacing  phantoms 

Spacing  phantoms  were  scanned  using  pCT  and  reconstructed  at  22  pm  isotropic  voxel  size.  These  images 
were  cropped  to  include  only  the  phantom,  and  were  resliced  to  expose  a  central  cross  section  representing 
the  phantom’s  banding  pattern.  DTS  scans  were  performed  on  the  phantom  tray  oriented  in  3  directions: 
parallel,  perpendicular,  and  at  22.5°  oblique  to  the  scanning  direction  (Figure  3).  As  there  were  2  sets  of 
spacing  phantoms  (oriented  in  line  and  at  22.5°  to  length  of  the  tray)  for  in-plane  assessment,  the  3  scanning 
directions  allowed  the  phantoms  to  be  scanned  a  total  of  4  different  angles  relative  to  the  scanning  direction 
(0°,  22.5°,  67.5°  and  90°  to  the  scanning  direction). 


(a)  (b)  (c) 

Figure  3:  DTS  reconstructions  of  spacing  and  thickness  phantoms  oriented  parallel  (a),  perpendicular  (b),  and 

22.5°  (c)  to  the  scanning  direction 
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Visual  Analysis  of  DTS  Spacing  Phantom  Images 

DTS  images  of  spacing  phantoms  were  visually  assessed  for  presence  of  a  banding  pattern.  A  central  cross- 
sectional  slice  from  each  phantom  was  viewed  at  200%  magnification  with  display  range  set  to  the  full  pixel 
value  range  of  the  image. 

Quantitative  Analysis  of  Microcomputed  Tomography  and  DTS  Spacing  Phantom  Images 
A  vertical  line  with  thickness  of  approximately  2  mm  was  drawn  between  the  top  aluminum  layer  and  the 
bottom  aluminum  layer  (Figure  4).  A  mild  median  filter  with  a  kernel  of  2  was  applied  to  the  image  to  reduce 
noise.  A  gray  value  profile  was  plotted  along  the  vertical  line,  which  was  processed  using  a  Fast  Fourier 
Transform  (FFT)  to  assess  frequency  of  the  banding  pattern. 

The  resulting  frequency-domain  plot  was  assessed  to  identify  maximum  peaks  representing  frequency 
magnitude  of  the  banding  pattern.  Average  layer  thickness  was  calculated  as  0.5*(1/frequency)  at  the 
maximum  peak.  The  peak  location  and  thickness  was  considered  a  gold  standard  by  which  to  compare 
subsequent  analyses  of  DTS  images.  The  same  procedure  was  performed  on  DTS  phantom  scans  from  a 
single,  randomly-chosen  scanning  session. 

Textural  information  was  visually  present  in  phantoms  250  pm  and  below  which  was  not  identified  in  FFT 
analysis,  so  mean  intercept  length  (MIL)  and  line  fraction  deviation  (LFD)  analyses  were  performed.  A  single 
cross-section  from  the  six  phantoms  at  the  0°  angle  was  cropped  from  the  DTS  reconstruction.  MIL  and  LFD 
analyses  were  performed  using  a  custom  MATLAB  script  (Sander  and  Barocas,  2009).  From  MIL  analysis, 
maximum  MIL  (MIL. Max)  and  its  angle  (MlL.MaxAngle)  were  recorded,  as  well  as  average  MIL  (MlL.Av), 
degree  of  MIL  anisotropy  (ratio  of  maximum  to  minimum  MIL,  MIL. DA),  and  standard  deviation  of  MIL 
(MIL.SD).  From  LFD  analysis,  maximum  LFD  (LFD. Max)  and  its  angle  (LFD.MaxAngle),  minimum  LFD 
(LFD. Min),  average  LFD  (LFD.Av),  degree  of  LFD  anisotropy  (ratio  of  maximum  to  minimum  LFD,  LFD. DA), 
and  standard  deviation  of  LFD  (LFD.SD)  were  recorded.  In  addition,  similar  parameters  were  recorded  when 
an  ellipse  was  fit  to  LFD  measurements  (LFD.EIlipseMax,  LFD.EIlipseMin,  and  LFD.EIlipseMaxAngle).  LFD 
and  MIL  parameters  allow  quantification  of  pattern  directionality  in  the  phantom  images.  LFD  is  typically  more 
sensitive  to  anisotropy  than  MIL,  as  the  method  calculates  the  fraction  of  bright  pixels  in  a  grayscale  image 
along  many  closely  packed  lines  at  various  angles,  rather  than  counting  intersections  at  various  angles  in  a 
binarized  image  (Geraets,  1998). 

Fractal  dimension  and  lacunarity  were  calculated  from  the  same  images  using  FracLac  (Karperien,  2007).  The 
box  counting  method  (Fazzalari  and  Parkinson,  1996;  Plotnick  et  al.,  1996)  was  used  to  assess  fractal 
dimension,  a  measure  of  image  complexity  which  uses  the  difference  between  maximum  and  minimum  gray 
values  to  describe  how  a  pattern’s  details  change  with  scale.  Lacunarity,  a  measure  of  translational  invariance, 
was  also  calculated.  Both  methods  work  together  to  characterize  the  existence  of  patterns  within  an  image. 
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Figure  4:  Creating  line  profile  for  Fast  Fourier  Transform  from  filtered  pCT  reconstruction  (upper  left).  In  the 
above  plot,  a  peak  was  identified  at  351 .23  pm  (upper  right).  The  same  procedure  was  repeated  for  DTS 

images  (bottom  center). 


Thickness  Phantoms 

Thickness  phantoms  were  scanned  using  microcomputed  tomography  and  reconstructed  at  10  pm  isotropic 
voxel  size.  These  images  were  cropped  to  include  only  the  phantom,  and  were  resliced  to  expose  a  central 
cross  section  representing  the  phantom’s  banding  pattern.  The  two  sets  of  in-plane  thickness  phantoms 
(oriented  in  line  and  at  22.5°  to  length  of  the  tray)  were  also  DTS  scanned  in  three  orientations  as  described 
above,  yielding  4  different  angles  relative  to  the  scanning  direction  (0°,  22.5°,  67.5°  and  90°;  Figure  3). 

Visual  Analysis  of  DTS  Thickness  Phantom  Images 

DTS  images  of  thickness  phantoms  were  visually  assessed  for  presence  of  aluminum  layers  separated  by 
plastic.  A  central  cross-sectional  slice  from  each  phantom  was  viewed  at  200%  magnification  with  display 
range  set  to  the  full  pixel  value  range  of  the  image. 

Quantitative  Analysis  of  Microcomputed  Tomography  and  DTS  Thickness  Phantom  Images 
The  same  FFT,  morphometric  and  textural  analyses  as  previously  described  were  performed  on  pCT  and  DTS 
images  of  the  thickness  phantoms. 


Aluminum  grid  panel  phantom 

The  aluminum  grid  panel  was  scanned  in  four  angular  orientations  relative  to  the  scanning  direction  (Figure  5; 
0°,  15°,  22.5°  and  45°).  Each  scan  generated  a  stack  of  1 1  reconstructed  image  slices  with  1  mm  spacing 
between  each  slice.  Two  rectangular  regions  of  interest  (ROI:  LI  and  L2)  with  500  pixel  width  and  maximum 
length  allowable  within  each  the  DTS  image  were  analyzed  for  the  gray  value  profiles,  covering  approximately 
60-70%  of  examined  image  area. 
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Figure  5:  Aluminum  grid  panel  was  imaged  at  four  angular  orientations  (0°,  15°,  22.5°  and  45°)  respective  to 
the  scanning  movement.  Each  angular  orientation  constituted  of  1 1  DTS  slices  that  were  spaced  1  mm 
between  slices.  The  gray  value  profile  of  each  image  slice  was  measured  at  two  principal  orientations  LI  and 

L2  (yellow  boxes). 


This  task  is  complete. 

Task  3.  Microcomputed  tomography  (pCT)  (Months  5-6): 

3a.  Microcomputed  tomography  scanning/reconstruction  of  individual  vertebral  bodies  at  30  micrometer  voxelsize 
using  a  calibrated  radiographic  reference.  (Month  5) 

3c.  Preparation  of  regions  of  interest,  segmentation  and  quantification  of  average  and  standard  deviation  of 
conventional  micro-CT  parameters  using  specialized  CT-analyzer  software:  bone  volume  fraction,  bone  surface 
to  volume  ratio,  trabecular  number,  thickness  and  separation,  structure  model  index,  connectivity  and  anisotropy. 
(Month  5) 

3d.  Preparation  of  regions  of  interest  and  calculation  of  fractal  dimensions  and  lacunarity.  (Months  5-6) 

3e.  Preparation  of  regions  of  interest,  segmentation  and  calculation  of  line  fraction  deviation.  (Month  6) 

We  have  previously  noted  that  calculations  of  microstructural  properties  are  affected  by  the  voxel  size  used  in 
pCT  scanning  and  image  reconstruction  (Kim  et  al.,  2004).  The  maximum  resolution  that  can  be  achieved  during 
a  pCT  scan  is  affected  by  the  size  of  the  object;  larger  vertebrae  allow  for  larger  minimum  voxel  size  (resulting  in 
lower  resolution).  In  order  to  avoid  a  potential  confounding  effect  of  vertebral  size  (through  voxel  size)  on 
calculated  parameters,  we  determined  a  minimum  voxel  size  that  can  be  used  for  all  vertebrae  in  the  study  by 
measuring  vertebral  dimensions  of  the  actual  specimens,  either  on  the  specimen  using  a  caliper  or  on  images 
captured  from  pilot  pCT  views.  We  determined  for  the  first  40  vertebrae  (T6,  T8,  Til  and  L3  vertebrae  from  5 
female  and  5  male  cadavers;  Age  63-90  years)  that  45  micrometers  is  the  minimum  voxelsize  that  can  be  used 
for  all  the  samples.  Large  scale  finite  element  models  of  cancellous  bone  based  on  pCT  images  with  50 
micrometer  voxelsize  achieved  good  agreement  with  experimentally  measured  stiffness  in  previous  work  (Hou  et 
al.,  1998)  indicating  that  the  cancellous  microstructure  of  vertebral  bone  is  sufficiently  resolved  at  this  resolution. 

The  in-house  pCT  system  used  in  the  project  is  built  on  the  principles  of  an  earlier  in-house  system  previously 
described  in  detail  (Reimann  et  al.,  1997).  The  presently  operating  system  uses  a  Kevex  16-watt  x-ray  source 
with  a  9-micron  focal  spot,  a  1888x1 408-pixel  Varian  PaxScan  2520  flat  panel  digital  x-ray  detector  with  127- 
micron  pitch,  a  Newport  precision  rotational  stage,  and  control  software  running  under  Windows  XP.  Each 
vertebra  was  placed  in  a  radiolucent  tube,  sealed  prior  to  scans,  and  scanned  using  our  in-house  pCT  system  at 
a  voxel  size  of  45  pm.  In  each  reconstruction  of  whole  vertebral  body,  the  largest  possible  cubical  volume  of 
interest  (VOI)  consisting  only  the  cancellous  bone  was  digitally  cut  out  using  ImageJ  (NIH,  USA).  The  bone  and 
non-bone  phases  were  segmented  using  Microview  (v2.0,  GE  Healthcare,  Little  Chalfont,  UK),  which  is  based  on 
standard  global  thresholding  procedures  previously  described  and  widely  used  (Kuhn  et  al.,  1990;  Hou  et  al., 
1998;  Yeni  and  Fyhrie,  2001). 

For  analysis  of  the  microstructure,  the  segmented  images  were  binarized.  The  analysis  was  carried  out  using 
both  custom-written  software  (Goulet  et  al.,  1994;  Yeni  et  al.,  2011)  and  commercial  software  (CT  Analyzer, 
Skyscan,  Kontich,  Belgium).  The  intra-specimen  mean  (Av)  and  standard  deviation  (SD)  of  the  following 
microstructural  properties  were  calculated  from  each  binarized  VOI  using  stereologic  methods:  Bone  volume 
fraction  (BV/TV),  bone  surface  to  volume  ratio  (BS/BV)  trabecular  thickness  (Tb.Th),  trabecular  number  (Tb.N) 
and  trabecular  separation  (Tb.Sp).  Degree  of  anisotropy  (DA)  was  also  calculated  using  the  3D  Secants, 
configuration  and  star  volume  methods.  The  Av  and  SD  were  recorded  as  measures  of  intraspecimen  magnitude 
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and  heterogeneity,  respectively,  of  a  given  parameter.  Connectivity,  structure  model  index  and  fractal  dimension 
parameters  were  calculated  using  CT  Analyzer  among  other  properties  properties  of  the  microstructure. 

We  have  automated  LFD,  MIL,  and  fractal  measurement  and  are  able  to  batch-process  image  sets  for  this 
analysis.  We  can  not  only  analyze  selected  sections  from  the  image  as  was  originally  planned  but  can  also 
process  an  entire  series  of  slices  within  a  single  pCT  volume.  Using  this  approach,  we  performed  LFD,  MIL  and 
fractal  analyses  using  volumes  of  pCT  images  matched  for  both  thickness  and  location  to  the  DTS  slice  (see 
Task  4). 

The  work  described  thus  far  fulfills  the  scope  of  the  statement  of  work. 

In  addition  to  the  analyses  of  microstructure,  we  performed  finite  element  simulations  on  each  of  the 
cancellous  bone  volumes,  in  order  to  estimate  cancellous  bone  stiffness  and  trabecular  stress  distribution 
properties.  This  was  to  obtain  Young’s  modulus  (EFem)>  coefficient  of  variation  of  von  Mises  stress  distribution 
(VMCV)  and  average  trabecular  von  Mises  stress  per  apparent  uniaxial  stress  (VMExp/aapp)  in  the  cancellous 
bone,  and  later  to  build  multi  regression  models  with  DTS  parameters  to  assess  prediction  capability  of  the 
DTS  system.  Linear  elastic  FE  models  with  isotropic  element  properties  were  prepared  using  a  custom-written 
meshing  algorithm  as  in  previous  studies  (Yeni  et  al.,  2005).  Young’s  modulus  of  each  bone  element  was 
assigned  using  an  approach  similar  to  that  used  by  Nekkanty  et  al.  (Nekkanty  et  al.,  2010)  resulting  in  the 
relationship: 

^tissue  (MPa)  =  (8020/1 130.257)  GLeiement 

where  GLeiement  is  the  gray-level  of  the  element  and  1130.257  is  the  gray  level  corresponding  to  a  mineral 
density  of  1.223g/cm3.  All  elements  were  assigned  a  Poisson’s  ratio  of  0.3.  Uniaxial  compression  was 
simulated  by  using  a  displacement  equivalent  to  an  apparent  strain  of  0.005.  The  models  were  displaced  in 
the  superoinferior  direction  while  the  displacements  in  other  directions  were  not  constrained  (Free-end 
conditions).  The  models  were  solved  using  a  custom-written  element-by-element  pre-conjugated  gradient 
solver  as  in  previous  studies  (Hou  et  al.,  1998;  Yeni  et  al.,  2005).  Apparent  modulus  of  cancellous  bone  was 
calculated  using  nodal  reaction  forces,  the  apparent  surface  area  of  the  sample  and  the  input  apparent  strain. 
The  average  and  standard  deviation  of  trabecular  von  Mises  stress  were  calculated  for  the  entire  sample. 

This  task  is  complete. 

Task  4.  Digital  Tomosynthesis  (DTS)  of  vertebrae  (Months  7-9): 

4a.  DTS  scanning  of  vertebral  bodies  using  the  SonialVision  Satire  II  system  (Shimadzu  Co.,  Kyoto,  Japan), 
with  the  superoinferior  axis  as  the  axis  of  image  rotation.  Scans  with  the  coronal  and  sagittal  views  centered  at 
zero  degree  angle  will  be  obtained.  (Month  7) 

4b.  Preparation  of  regions  of  interest  and  calculation  of  fractal  dimensions  and  lacunarity.  (Months  7-8) 

4c.  Preparation  of  regions  of  interest,  segmentation  and  calculation  of  line  fraction  deviation.  (Month  8) 

Using  the  Shimadzu  Sonialvision  Safire  II  tomosynthesis  system  (Shimadzu,  Kyoto),  74  radiographic 
tomosynthesis  images  were  captured  over  2.5  seconds  at  65  kV  and  250A  and  then  reconstructed  with  image 
slice  thickness  of  1  mm. 

We  conducted  our  study  in  two  phases  to  fully  explore  the  potential  of  DTS  in  predicting  fracture  strength  of 
vertebral  bones.  The  first  phase  was  conducted  as  a  preliminary  study  to  test  the  feasibility  of  MIL,  LFD  and 
fractal  analyses  on  DTS  images  and  to  identify  which  parameters  from  these  analyses  would  be  informative  in 
predicting  cancellous  bone  stiffness  obtained  from  pCT-based  finite  element  analysis.  In  phase  1,  imaging 
parameters  were  established  and  image  pre-processing  techniques  were  examined. 

In  the  phase  2  of  the  study,  with  an  enhanced  understanding  of  which  DTS  parameters  are  informative  in 
stiffness  prediction,  we  then  investigated  the  effect  of  scanning  orientation  and  3D  analysis  of  DTS  slices 
predicting  both  pCT  stereological  parameters  and  mechanical  fracture  strength  of  the  40  vertebral  bones. 
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Phase  1 

In  phase  1,  T6,  T8,  Til,  and  L3  vertebral  bodies  were  scanned  in  thin-walled  plastic  cylinders,  aligned  and 
separated  using  radiolucent  foam.  The  cylinders  were  scanned  transverse  to  the  axis  of  the  spine  with  the 
cylinders  aligned  such  that  reconstructed  slices  represented  the  anterior-posterior  (AP)  anatomic  direction.  A 
single  coronal  image  from  each  bone  was  selected  for  analysis.  The  images  were  then  processed  with  three 
different  image  processing  methods  (Figure  6)  in  an  effort  to  enhance  extraction  of  trabecular  structures  within 
the  image  texture.  Three  types  of  processing  methods  were  examined: 

i)  USM  -  unsharp  mask  (USM)  filter  sharpened  without  binarization  (i.e.,  grayscale) 

ii)  Bin.Csv  -  conservative  binarization  by  use  of  weak  convolution  and  high  pass  filters 

iii)  Bin.Ag  -  aggressive  binarization  by  use  of  strong  convolution  and  high  pass  filters 

The  filtering  processes  were  performed  using  Photoshop  (Adobe  Systems  Inc,  San  Jose,  California)  and  ImageJ 
(NIH,  Bethesda,  Maryland).  However  fractal  analysis  was  performed  on  images  with  no  pre-processing. 


Figure  6:  A)  A  typical  DTS  image  segment  used  for  the  analyses.  There  were  two  types  of  ROI;  B)  Fractal 
dimension,  lacunarity  and  S^  were  calculated  through  fractal  analysis  without  image  pre-processing.  C)  MIL 
and  LFD  parameters  (Max,  DA,  Av  and  SD  of  MIL  and  LFD)  were  calculated  for  each  of  the  three  types  of 

image  pre-processing  methods  (USM,  Bin.Csv  and  Bin.Ag). 

Fractal  analysis  -  Using  box  count  fractal  analysis  method  (Fazzalari  and  Parkinson,  1996;  Plotnick  et  al., 
1996),  fractal  dimension  (FD)  and  mean  lacunarity  (X)were  obtained  by  use  of  FracLac  2.5  plugin  (Karperien, 
2007)  of  ImageJ  (NIH,  Maryland).  The  box  count  method  uses  a  series  of  sequentially  reducing  box  size  grids 
on  the  image  to  sample  the  local  difference  between  the  maximum  and  minimum  gray  values  in  each  box  region. 
As  a  secondary  derivative  measure  of  lacunarity,  a  slope  of  lacunarity  against  box  size  value  was  also  calculated 
(S„). 

MIL  and  LFD  analyses  -  We  have  performed  both  MIL  and  LFD  analyses  using  custom  MATLAB  codes 
(Mathworks,  MA),  calculating  MIL  and  LFD  values  for  each  1  degree  increment  over  360  degrees  (Sander  and 
Barocas,  2009).  The  MIL  and  LFD  values  were  polar  plotted  (Figure  7)  and  mean  (Av),  standard  deviation  (SD), 
maximum  (Max)  and  degree  of  anisotropy  (DA)  of  MIL  and  LFD  were  obtained  for  each  DTS  image. 


LFD  (line  fraction  deviation).  Note  how  overall  orientations  of  the  polar  plots  correspond  to  that  of  the 
trabecular  texture  in  (A).  Dotted  red  lines  represent  the  principal  components  in  MIL  and  LFD.  The  degree  of 
anisotropy  (DA)  was  calculated  from  the  ratio  of  maximum  to  secondary  principal  component  where  MIL. DA 
was  from  MlL.Max/MlL.Min  and  LFD. DA  was  from  LFD. Max/LFD. Ortho. 


In  order  to  further  understand  LFD,  MIL  and  fractal  parameters,  we  subjected  a  single  DTS  slice  (1  mm  thick) 
and  20  corresponding  pCT  slices  (each  0.045  mm  thick)  -  the  regions  matched  for  both  thickness  and  location 
-  to  LFD,  MIL  and  fractal  analyses  (Figure  8)  and  constructed  pair-wise  mixed  regressions. 


Figure  8:  DTS  slice  (a)  and  pCT  stack  of  equivalent  thickness  (b)  were  analyzed  for  MIL,  LFD  and  fractal 
parameters  and  then  DTS  and  pCT  generated  parameters  were  compared.  DTS  image  has  been  inverted  for 
white  bone  presentation. 


The  work  described  thus  far  fulfills  the  scope  of  the  statement  of  work. 

During  this  phase,  we  developed  the  capability  to  automatically  batch-process  a  large  number  of  images  which 
allows  us  to  analyze  multiple  slices  of  a  single  DTS  reconstruction.  In  analyzing  the  first  set  of  images,  it  became 
apparent  that  the  resolution  of  trabecular  structure  is  highly  directionally  dependent  in  DTS  imaging.  It  was  also 
found  that  image  pre-processing  does  not  conclusively  enhance  correlation  between  pCT  stereology  parameters 
and  DTS  texture  parameters  and  was  therefore  not  used  in  future  studies.  With  the  knowledge  gained  from 
phase  1,  we  were  posed  to  better  address  the  research  topic. 

Phase  2 

With  the  better  understanding  of  which  DTS  parameters  are  potentially  informative  in  predicting  cancellous  bone 
finite  element  stiffness  from  phase  1  of  the  study,  we  rescanned  the  same  40  vertebral  bones  but  this  time  in  six 
different  configurations  (Figure  10);  in  coronal  (AP  -  anterior  posterior)  and  sagittal  (LM  -  lateral  medial)  views  at 
three  different  angles  to  the  spinal  axis  (par,  ang,  perp  -  0°,  22.5°  and  90°)  so  as  to  understand  which  view  and 
orientation  would  provide  best  prediction  in  fracture  strength. 

Precision  in  specimen  positioning  was  considered  an  issue  in  phase  1,  so  specimens  were  mounted  and 
consistently  aligned  in  a  custom  radiolucent  scanning  tank  filled  with  0.9%  saline  (Figure  11).  Alignment  was 
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ensured  using  a  radiolucent  clamping  system. 


Figure  10:  Vertebrae  were  scanned  at  six  orientations  in  total,  parallel  (a),  oblique  (b),  and  perpendicular  (c) 
to  the  superoinferior  axis  of  the  spine,  in  both  AP  (anterior-posterior)  and  LM  (lateral-medial)  views. 
Superior/inferior  (S/I)  orientation  was  aligned  with  respect  to  the  scanning  direction  (indicated  by  large  arrows). 


Figure  11:  Specimens  were  mounted  on  trays  inside  custom-built  radiolucent  tanks, 
allowing  for  precise  specimen  positioning  and  orientation  with  respect  to  image  axes  (a). 

Consistent  AP/LM  alignment  and  transverse  tilt  correction  was  achieved  via  rotation  on  a 
tray-mounted  clamping  system  within  each  box  (b). 

Matching  cuboidal  VOIs  were  prepared  from  both  pCT  and  DTS  reconstructions.  The  DTS  volume  was  analyzed 
on  a  slice-by-slice  basis  for  texture  parameters  (LDF,  MIL,  and  fractal)  and  the  pCT  volume  was  analyzed  for 
stereology  parameters.  A  description  of  the  measured  parameters  and  nomenclature  are  described  in  Figure  12. 
Mixed  regressions  were  used  to  compare  DTS  (LFD,  MIL  and  fractal)  parameters  to  pCT  stereology  parameters. 
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Nomenclature 

X(Y.Z) 

^  1^ 

Between  planes:  Analysis: 
Av,  average  MIL 

SD,  std  dev  LFD 

F,  Fractal 


Within  plane 

parameter: 

Av 

SD 

Max,  maximum 
DA,  anisotropy 
D,  dimension 
A,  lacunarity 
S/„  lacunaritv/scale 


Figure  12:  (A)  A  square  region  cropped  out  of  the  DTS  image  is  analyzed  using  the  fractal,  MIL  and  LFD 
analysis  methods.  MIL  and  LFD  are  calculated  in  all  directions  in  the  image.  (B)  A  typical  polar  plot 
representation  of  the  LFD  data,  from  which  an  average,  standard  deviation,  max  and  anisotropy  (LFD  at  90°  / 
LFD  at  0°)  can  be  calculated  for  each  slice.  (C)  Each  slice  in  the  stack  is  analyzed.  Stack  average  and 
standard  deviation  (e.g.  SD  (LFD.SD))  are  calculated  as  measures  of  3D  microstructure  within  the  bone. 


Additionally  we  have  measured  approximate  vertebral  height  and  width  using  the  central  slices  from  AP- 
transverse  and  LM-transverse  (corresponding  to  coronal  and  sagittal  views  -  Figure  13) 
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Figure  13:  A  central  slice  from  (a)  anterior-posterior  (AP)  coronal  and  (b)  lateral-medial  (LM)  sagittal 
scans  was  used  to  measure  waist-to-waist  width  and  endplate-to-endplate  height. 


In  phase  2,  we  also  investigated  three  methods  of  ROI/VOI  selection  (Figure  14)  to  determine  if  increasing 
analysis  volume  would  improve  the  prediction  performance.  The  three  methods  were  (A)  using  a  single  largest 
possible  ROI  from  AP  and  LM  mid-plane  (a  single  2D  image),  (B)  multiple  slice  VOI  with  constant  dimension 
square  within  each  slice  and  (C)  multiple  slice  VOI  with  largest  possible  square  ROI  in  each  given  slice 
resulting  in  a  varying  cross-sectional  ROI  area  across  the  stack.  At  the  end  of  phase  2,  we  constructed  both 
mixed  regression  models  predicting  peak  fracture  forces  of  the  40  vertebrae  from  wedge  compression  tests. 
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Analysis  area/volume 


AP  DTS  slice 
LM  DTS  slice 


Figure  14:  Three  ROI  methods  in  selecting  DTS  images  for  analysis.  (A)  single  largest  possible  ROI  from  AP 
and  LM  mid-plane,  (B)  VOI  with  constant  dimension  or  (C)  VOI  with  largest  square  ROI  per  slice. 


This  task  is  complete. 


Task  5.  Mechanical  testing  of  vertebral  bodies  (Months  10-11) 

5a.  Non-uniform  compression  testing  of  vertebral  bodies  using  a  materials  test  machine  (8501 M,  Instron  Inc., 
MA).  Determine  vertebral  strength  as  the  maximum  load  during  fracture.  (Months  10-11) 

5b.  Reduction  of  mechanical  test  data.  (Month  11) 

Preliminary  wedge  loader  model 

The  mechanical  testing  fixture  and  protocol  were  designed  to  apply  non-uniform  loading  conditions  to  the 
specimen  vertebral  bodies.  Specifically,  the  non-uniform  loading  would  consist  of  uniaxial  compression  and  a 
bending  moment,  which  would  simulate  the  loads  producing  a  wedge  fracture,  often  a  result  of  heavy  weight 
lifting.  The  objective  of  the  design  and  protocol  was  to  allow  independent  control  of  the  compressive  load  and 
the  bending  moment.  Previous  researchers  achieved  combined  loading  via  wedged  load  platens  (McDonnell  et 
al.,  2010),  multi-arm  hydraulic  systems  (Whealan  et  al.,  2000),  and  robotic  arm  (Buckley  et  al.,  2007).  We 
produce  the  requisite  loading  via  an  axial  and  rotary  actuator.  The  axial  actuator  control  applies  the 
compressive  loading  condition  while  the  rotary  actuator  control  applies  the  bending  loading. 

During  the  study  period,  we  started  by  designing  and  fabricating  a  prototype  combined  loading  fixture  to  prove 
out  the  concept  -  Figure  15  is  a  conceptual  design  rendering  of  the  prototype.  The  upper  load  applicator  has 
an  inclined  surface  which,  when  rotated,  directs  one  of  the  load  transfer  rods  downward  and  the  other  upward. 
This  rotates  the  combined  loading  platen  and  imparts  bending  to  the  specimen  mounted  to  the  specimen 
mount.  Axial  displacement  of  the  load  applicator  directs  both  load  transfer  rods  downward  simultaneously  and 
imparts  specimen  compression.  Testing  during  the  prototype  trials  consisted  of  applying  compression, 
bending,  and  combined  compression  and  bending  to  a  foam  test  block  to  evaluate  the  combined  loading  fixture 
operation.  The  combined  loading  of  the  test  foam  blocks  produced  a  linear  load  response  (Figures  16-19). 
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Figure  15:  Conceptual  design  rendering  of  the  combined  loading  fixture 

Fabrication  of  the  project  combined  loading  fixture  commenced  after  the  prototype  tests  were  complete 
(Figures  20-22).  Two  design  changes  were  incorporated  in  the  modified  design.  The  first  design  change 
consisted  of  reducing  the  angle  of  the  inclined  surface  on  the  load  applicator  from  22°  to  1 1 0  which  reduced  the 
maximum  bending  angle  placed  on  the  vertebral  body  to  5°.  The  second  design  change  was  to  include  a  flat 
section  on  half  of  the  load  applicator,  with  the  other  half  inclined  at  the  desired  angle  noted  above.  The 
purpose  of  the  flat  section  on  the  load  applicator  was  to  prevent  lift  off  of  the  combined  loading  platen  from  the 
posterior  rim  of  the  vertebral  body. 
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Foam  Test  Block  -  Test  2 
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Figure  16:  Prototype  test  #2.  Foam  block  spinal  surrogate. 
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Figure  17:  Prototype  test  #3.  Foam  block  spinal  surrogate. 
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Foam  Test  Block  -  Test  4 
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- Load  -  Compression  Actuator  Rotation  - Actuator  Displacement 


Figure  18:  Prototype  test  #4.  Foam  block  spinal  surrogate. 


Foam  Test  Block  -  Test  5 
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—  Load  -  Compression  —  Actuator  Rotation  —  Actuator  Displacement 
Figure  19:  Prototype  test  #5.  Foam  block  spinal  surrogate. 
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Figure  20:  Preliminary  wedge  loading  fixture.  Load  applicator  attached  to 
the  test  machine  actuator  is  shown  with  flat  and  inclined  sections. 
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Figure  21 :  Preliminary  wedge  loading  fixture  -  lower  frame.  Lower  frame  shows  the  load 
transfer  rods,  loading  platen  and  specimen  mounted  and  attached  to  load  cell. 
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Figure  22:  Preliminary  wedge  loading  fixture  -  lower  frame,  loading  platen  and  specimen  mount.  Picture 
shows  a  Sawbones®  surrogate  vertebral  body  specimen  positioned  in  the  loading  fixture  between  the  loading 
platen  and  specimen  mount.  View  is  in  the  A-P  direction  and  bending  will  be  applied  anteriorly  in  the  sagittal 
plane. 


In  our  calculations  we  have  assumed  that  the  vertebral  body  loading,  with  the  combined  loading  fixture,  will 
conform  to  a  beam  theory  cantilever  model  (Figure  23).  The  model  assumes  a  fixed  end  (the  potted  inferior 
vertebral  endplate)  and  a  guided  opposite  end  (superior  vertebral  endplate).  The  input  forces  are  P 
(compressive  load)  and  W  (transverse  load),  with  W  placed  at  a  distance  a  from  the  superior  endplate.  In  this 

case,  a  equals  zero.  The  reaction  moments  at  the  superior  and  inferior 
P  endplates  can  be  calculated  from  these  loads.  The  forces  (P  &  W)  are 

directly  measured  with  a  3-axis  load  cell  (FUTEK  Advanced  Sensor 
Technology  Inc.,  Irvine  CA,  USA). 

Multiple  tests  were  conducted  on  a  Sawbones®  vertebral  body  using  the 
combined  loading  fixture.  Test  modes  included  compression  only,  bending, 
and  combined  compression  and  bending.  Compression  and  combined 
loading  tests  were  conducted  with  an  endpoint  uniaxial  displacement  of 
2mm  and  1mm,  respectively,  in  the  vertical  (z-axis)  direction  at  a  rate  of  0.1 
mm/s.  Bending  and  combined  loading  tests  were  conducted  with  a  rotary 
actuator  endpoint  rotation  of  30°-40°  and  20°-30°,  respectively,  at  a  rate  of 
1°/s.  The  observed  loads  were  within  the  expected  range  of  loads  for 
cadaveric  vertebrae.  Loads  ranged  from  4993N  to  9047N.  During  these 
tests  it  was  noted  that  the  z-axis  displacement  did  not  travel  to  the 
programmed  endpoints  or  at  the  programmed  displacement  rate.  The  result 
was  endpoint  displacements  that  increasingly  deviated  from  the 
programmed  endpoint  as  the  distance  travelled  increased  (Figure  24).  Also 
Figure  23:  The  model  for  the  the  displacement  behavior  at  the  transition  between  downward  and  upward 

multiaxial  (wedge)  loading.  travel  became  nonlinear.  The  first  three  displacement-time  curves  (1 .0  mm 

endpoint  -  Imm/s,  1.5  mm  endpoint  -  Imm/s,  2.0  mm  endpoint  -  Imm/s)  show  that  the  actual  displacement 
was  83%  and  77.0%  of  the  programmed  endpoint  for  1.0  mm  of  travel  and  2.0  mm  of  travel,  respectively.  The 
displacement  rates  had  similar  deviation  from  the  programmed  value.  This  issue  was  found  to  be  a  result  of  the 
PID  gain  settings  of  the  displacement  control  feedback  loop  in  the  materials  test  machine.  The  PID  gain 
settings  were  modified  to  optimum  settings  for  the  test  setup,  which  corrected  the  endpoint  displacement  and 
rate  to  a  maximum  of  91%  of  the  programmed  values  (1.0  mm  endpoint  -  0. Imm/s  -  P30I20),  without  the 
actuator  control  becoming  unstable. 
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Figure  24:  Comparison  of  z-axis  displacement  deviation  from  programmed  endpoint  and  displacement  rate. 
The  deviation  between  the  programmed  and  actual  travel  increases  as  a  function  of  the  endpoint  distance. 

Another  aspect  of  the  displacement  deviation  that  has  to  be  accounted  for  is  the  machine  compliance.  During 
a  mechanical  test,  the  recorded  displacement  of  the  actuator  is  comprised  of  specimen  and  test  apparatus 
deformations  (Kalidindi  et  al.,  1997): 


Sr  -  §c  +  Ss 

where,  §R  is  the  recorded  displacement  of  the  actuator,  8C  is  the  deformation  of  the  test  apparatus,  the 
combined  loading  fixture  in  this  case,  and  8S  is  the  specimen  deformation.  Tests  were  conducted  to 
characterize  the  machine  compliance  by  defining  an  expression  for  the  test  apparatus  deformations  in  terms  of 
the  applied  force  8C(F).  Upon  inspection  of  applied  displacement  axial  displacement  and  rotation  vs.  linear 
variable  displacement  transducer  (LVDT)  measurement  at  the  superior  platen  surface,  it  was  determined  that 
excessive  compliance  was  present  in  this  construct  due  to  unavoidable  internal  movement  within  the  linear 
bearing  and  ball  transfer  housing.  Using  the  information  gained  from  the  preliminary  wedge  loader  model,  a 
simplified  construct  was  fabricated  which  allowed  combined  loading  via  a  lever  arm  with  an  offset  hinge. 

Final  wedge  loader  model 

The  final  wedge  loader  model  produced  uniaxial  compression  and  a  bending  moment  through  the  use  of  the 
linear  actuator  only  (Figure  25).  In  this  construct,  a  specimen-specific  total  displacement  was  applied  at  the  top 
surface  of  a  flat  beam  lever  at  a  specimen-specific  distance  from  a  fixed  hinge  via  a  ball  transfer.  The  entire 
construct  is  free  to  translate  left  and  right  in  order  to  set  a  hinge-specimen  distance,  which  can  be  locked  to  the 
table  via  four  bolts.  The  height  of  the  hinge  point  is  set  before  each  test  such  that  the  loader  beam  is  level. 
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Figure  25:  Final  wedge  loader.  The  simplified  design  allows  for  accurate 
specimen  positioning  with  minimized  system  compliance. 


Figure  26:  Close-up  view  of  deformed  human  vertebral  body  in  final  wedge  loader  system. 


Forces  are  measured  at  the  inferior  end  via  a  FUTEK  triaxial  load  cell,  as  well  as  the  Instron  uniaxial  load  cell 
at  100  Hz  recording  frequency.  High  quality  video  was  also  recorded  throughout  testing  (Figure  26).  The  total 
displacement  was  calculated  such  that  the  test  introduced  a  mild  wedge  fracture  with  10%  difference  in  height 
between  posterior  and  anterior  heights  and  an  overall  25%  mean  reduction  in  vertebral  height  (see  calculations 
below).  Displacement  rate  was  constant  at  0.1  mm/sec. 
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Final  wedge  loader  calculations 


H 


W/2 - -  W/2 


Figure  27:  Combined  loading  parameters  were  defined  by  distance  between  the  specimen  center 
and  the  hinge  point  (R),  as  well  as  total  actuator  displacement  ( d ). 

In  order  to  introduce  a  pre-determined  ratio  of  anterior  height  (Ha)  to  posterior  height  ( Hp )  {HaIHp= 0.9)  and  a 
mean  specimen  reduction  in  height  c=0.75,  distance  between  specimen  center  and  hinge  point  ( R )  and  total 
actuator  displacement  ( d)  were  determined  using  the  undeformed  heights  (H)  and  anterior-posterior  widths  (W) 
of  each  specimen,  which  were  determined  from  CT  images.  The  specimen  was  potted  on  both  inferior  and 
superior  surfaces,  and  inferior  surface  was  fixed  while  the  superior  surface  was  free  to  translate  along  the 
loader  beam.  The  system  was  thus  modeled  with  a  polar  approach,  where  superior  contact  surfaces  were 
considered  to  translate  and  rotate  together  with  the  loader  beam  along  a  radius  originating  at  the  hinge  point. 

From  the  geometry  in  Figure  27,  the  following  constraints  were  established: 

(1)  L  =  R  +  A 

(2)  Lsind  =  d 

(3)  L  cos  9  =  R 

Substituting  L  from  (1)  into  equations  (2)  and  (3)  yielded: 

(4)  ( R  +  A)  sin  9  =  d 

(5)  (R  +  A)  cos  9  =  R 

Pythagorean  Theorem  was  used  to  find  an  expression  for  A  in  terms  of  R  and  d: 

L2  =  (R+  A)2  =  R2  +  d2 
R2  +  2RA  +  A2  =  R2  +  d2 
A2  +  2RA  -  d2  =  0 


The  y  component  of  A  was  calculated  as  such: 


(7)  hmy  =  Asin  9  =  Hm  —  (H  —  d) 

Substituting  the  positive  solution  for  (6)  for  A  and  equation  (4)  for  sin  0  yielded  the  following  relationship: 

(8)  (-R  +  VWTd?)  (^p)  =  Hm-(H-d) 
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Equation  (8)  was  solved  for  d: 

1  J  V(R2_(W_Hm)2) 


Hm  was  assumed  to  be  a  constant  times  the  undeformed  height: 
(10) 


Hm  c  x  H 


Equation  (10)  was  substituted  for  Hm  in  equation  (9): 

j  _  R(H-cH )  _  RH (1— c) 


d  = 


V(R2-(H-cH)2) 


J(r2-(h(  i-c))2) 


Considering  that  the  lines  drawn  by  distances  dp  and  da  form  similar  triangles  to  the  triangle  drawn  with 
distance  d,  (H-Hp)  and  ( H-Ha )  replaced  ( H-Hm )  in  equation  9  in  the  expressions  for  dp  and  da,  respectively. 
Likewise,  (R-W/2)  and  (R+W/2)  replaced  R  in  equation  9  in  the  expressions  for  dp  and  da,  respectively: 


(12) 

d  —  - 

(R-f  )(H-Hp) 

r  J 

(0 r-t)2-("-«p )2) 

(13) 

da  =  ~ i 

(R+y)(H-Ha) 

From  the  geometry  in  Figure  27,  we  know  the  following  relationships  hold: 
(14) 


—  —  =  tan  6 


R~r  ~ 

2  2 


Substituting  equations  (12)  and  (13)  for  dp  and  da  in  equation  (14)  yielded: 

("~"p)  _  in-ng) 


It  was  desired  that  the  mechanical  test  continued  until  ratio  Hp/Ha  reached  a  physiologically  relevant  constant  k\ 


(16)  Hf=k  ^  Hp  =  kHa 


Hp  from  equation  (16)  was  substituted  into  equation  (15)  and  then  solved  for  Ha\ 


(17) 


Ha  = 


HW 


HW 


kR+~R+^j  R(k-l)+y(k+l) 


The  similar  triangles  formed  by  distances  H-Ha  and  H-Hm  were  used  to  calculate  R: 
(18)  sin  9  =  -  JIzHsl- 

v  J  R  R+W/2 

Substituting  Hm=cH  into  equation  (18)  and  solving  for  R  yielded: 

C»>  «  =  (t)  (“) 

The  value  for  R  found  in  equation  (19)  was  substituted  into  equation  (17): 

(20)  (k  +  1  )Ha2  -  (k  +  1  +  2 d)HHa  +  2 cH2  =  0 

Solving  for  Ha  yielded  the  following  relationship: 

(21) 


j,  H  (k+l+2c)±V  (k+l+2c)2-8c(k+l) 

nn  =  —  X - 


(k+l) 

The  radicand  in  equation  (21)  must  be  equal  to  or  greater  than  0: 
(fc  +  l)2  +  4 c(k  +  1)  +  4c2  —  8 c(k  +  1)  >  0 
(k  +  l)2  -  4 c(k  +  1)  +  4c2  >  0 
[(k  +  1)  -  2c]2  >  0 
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Thus, 

( k  T  1)  —  2 c  ^  0  and  ( k  T  1)  —  2 c  <  0 

k>2c  —  l  and  k  <  2c  -  1,  therefore  all  values  of  k  and  c  are  admissible. 

The  value  found  for  Ha  in  equation  (21)  was  substituted  into  equation  (19),  forming  the  expression  for  R: 


(22) 


R 


=  (t)I 


H-cH 


H  (k+l+2c)±V (k+i+2c)2-8c(k+l)^ 
2X  (fe+1) 


cH- 


The  value  found  for  R  in  equation  (22)  was  substituted  into  equation  (11),  forming  the  expression  for  d: 


m 


H-cH 


(23) 


d  = 


„  ,  H  (k+l+2c)±V(k+l+2c)2-8c(fc+l)'\ 

CH~[  2X - (k+1) - J 


H(l-c) 


(( 


m 


i\\ 


Hw(k+i+2c)+V(k+l+2c)2-8c(k+l) 


(fc+l) 


-(H(l-c))2 


J 


Simplifying  this  expression  yielded: 
(24)  d  = - 


(c—l)2H(k+l)W 


((2c-l)k+V(-2c+fe+l)2-l) 


(c- 1)2 


y((2c-l)k+J  (- 2c+fc+l)2-l) 


In  addition,  /cand  cmust  be  chosen  appropriately  to  yield  Rand  d values  that  simultaneously  fulfill  several 
physical  constraints  of  the  mechanical  testing  setup  (Figure  28): 


Figure  28:  Physical  constraints  of  the  mechanical  testing  system  imposed  by  loader  geometry. 
Values  for  R  and  d  must  be  chosen  according  to  these  constraints  so  that  the  test  is  feasible. 


Constraint  1 :  If  the  load  cell  stack  is  moved  as  far  left  as  possible,  the  value  for  R  is  inadmissible  if  less  than 
the  radius  of  the  load  cell  stack  plus  half  the  thickness  of  the  vertical  support  which  contains  the  hinge  pivot 
point  at  its  center: 


R  > 


•>(¥H 

/126.85  mm\ 


W, 


HINGE 


- 


+ 


2  ) 
'19.1  mm\ 


R  >  72.975  mm 
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Constraint  2:  Taking  into  account  translation  of  the  loading  point  down  the  length  of  the  beam  during  rotation, 
the  value  for  R  must  be  less  than  the  original  length  of  the  beam  less  the  x  component  resulting  from  rotation 
minus  half  the  width  of  the  specimen: 

Using  known  quantities,  solve  for  theta: 


sin  9  = 


H- Hr, 


R 

6  =Sm  \  R  / 


Substitute  the  value  for  theta  found  above  into  the  following  expression  and  solve  for  WBeam’: 

„  Wbeam' 
cos  9  =  — - 

WBEAM 

Wbeam  =  WBeam  cos  9  =  195.5  mm  x  cos  9 

As  discussed  above,  the  following  relationship  must  hold: 

,  (W\ 

R  <  Wbeam'  ~  (y) 

(25)  R  <  Wbeam  cos  9 -(^f) 


W 

R  <  195.5  mm  x  cos  9 - 

2 

Constraint  3:  So  as  to  avoid  damaging  the  upper  right  corner  of  the  load  cell  with  the  free  end  of  the  loader 
beam,  the  vertical  displacement  of  the  loading  beam  at  a  distance  R+WLC/2  (called  DRCorner)  must  be  less 
than  the  undeformed  specimen  height. 

From  the  geometry  drawn  above,  we  know  the  following  relationships: 

d 

ta„9=- 

DrcORNER 


tan#  = 


W, 


LC 


R  + 

£ 

Set  these  expressions  equal  to  each  other  and  solve  for  DRCorner: 

p  _  Prcorner 
R  ~ " 


R  + 


WLC 


2 

cl  /  IFlc  ' 
Hrcorner  =  —  X  ( R  +■ 


h'f)- 


R  \  2 

As  discussed  above,  the  following  relationship  must  hold: 

(26) 


d  +  ■ 


dW, 


LC 


2  R 


d  +  d-^<H 

2R 


Constraint  4:  The  pot  has  a  3”  inside  diameter  (DPOt)-  The  Feret  diameter  (DF)  was  measured  from  a  max 
intensity  projection  of  all  axial  CT  slices  for  a  binarized  mask  of  each  vertebra.  The  Feret  diameter  is 
essentially  a  “caliper  diameter”  (i.e.,  the  diameter  of  the  smallest  circle  that  would  completely  enclose  the 
shape). 

Dp  <  Dpot 

This  constrain  held  true  for  the  largest  vertebra,  thus  the  3”  pot  was  appropriate  for  all  specimens. 


Constraint  5:  To  facilitate  the  hinge  mechanism,  the  loader  beam  has  a  slot  milled  in  it  which  extends  1 .91  cm 
from  the  vertical  support.  Through  rotation,  the  posterior  edge  of  the  superior  construct  surface  (a  distance  RP 
from  the  pivot  point)  cannot  slide  onto  the  slotted  portion  of  the  loader  beam,  so  RP  in  the  deformed 
configuration  must  be  more  than  half  the  hinge  width  (WHinge/2)  plus  the  slot  length  (WSlot)-  This  constraint  can 
be  expressed  in  terms  of  R\ 

Where  RP  =  R 


CL  and  RP  >  WSL0T  +  CLsmss. 


n  W  ^  w  ,  W 'hINGE 

R  ~  y  >  WSL0T  H  - 
(27)  R>Wslot+^ f^  +  7 
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Equations  (22)  and  (24)  were  solved  for  Ft  and  d,  respectively,  for  an  array  of  possible  k  and  c  combinations  (in 
0.025  increments)  using  specimen-specific  geometry  values  W and  /-/determined  from  CT  images.  This  was 
done  for  the  smallest  specimen  (lowest  1/1/ and  H)  as  a  worst  case  scenario,  because  it  would  produce  the 
lowest  /?  value.  The  combination  of  /rand  c was  simultaneously  optimized  such  that  a  physiologically  relevant 
wedge  fracture  was  produced  while  system  constraints  1-5  were  fulfilled.  Optimized  /cand  c  values  were 
established  via  this  method  and  applied  for  all  specimens,  such  that  the  mechanical  test  would  end  at  25% 
mean  reduction  in  height  with  10%  difference  between  posterior  and  anterior  heights.  Subsequently  solving 
equations  (22)  and  (24)  to  establish  R  and  d  values  for  all  40  specimens  fulfilled  constraints  1-5  for  all 
specimens. 

The  40  human  vertebral  bodies  were  tested  to  failure  following  all  imaging  studies.  An  example  failure  curve  is 
shown  in  Figure  29.  Local  maximum  at  the  first  peak  was  considered  failure  load. 


Combined  Loading  of  Human  L3  Vertebra 


Figure  29:  Example  force-displacement  curve  from  a  human  L3  vertebra  using  the  final  wedge  loader  system. 
Failure  load  was  recorded  for  the  purpose  of  correlation  with  texture  parameters  from  tomosynthesis  images. 
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Results: 


Results  from  Task  2 


Density  Phantom  Images 

If  the  analysis  region  is  limited  to  the  central  portion  of  the  cylinder,  a  nonlinear  relationship  was  found  between 
radiographic  density  (DTS  gray  values)  and  solution  density  (Figure  30a,  b)  (A  cubic  polynomial  explains  over 
99.9%  of  this  relationship).  The  within-session  and  between  session  variability  of  density  phantoms  was  less 
than  5%  (coefficient  of  variation=0.5%  to  4.7%)  with  largest  variations  caused  by  air  (3.5%  and  4.7%  for  within- 
and  between-session  variability,  respectively). 


Radiograph  gray  value  was  found  to  have  a  linear  response  to  increased  solution  density  (Figure  31). 
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Figure  30:  (a)  Average  DTS  gray  value  versus  solution  density  for  5  sessions,  (b)  Average  DTS  gray  value 
versus  solution  density  within  a  single  session.  Error  bars  represent  between-session  gray  value  range. 
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Figure  31 :  Average  DTS  gray  value  versus  solution  density  vs.  solution  concentration  obtained  from  DTS 
radiogram.  Error  bars  represent  the  range  of  gray  values.  (R2=0.9981  for  linear  regression.) 

Although  the  nonlinear  relationship  found  between  gray  values  and  density  in  DTS  images  may  be  useful  for 
quantitative  purposes,  the  edge  enhancement  apparent  in  DTS  images  may  affect  neighboring  radiographic 
density  values  making  quantitative  density  measurements  difficult.  On  the  other  hand,  the  high  resolution 
digital  radiography  capabilities  of  the  DTS  system  offer  another  possibility  for  density  measurements  (Kinds  et 
al.,  2011).  Combined  with  the  structural  analyses  from  the  DTS  images,  density  measurements  form  the  DTS 
radiograms  may  result  in  a  powerful  tool  for  bone  quality  assessment. 

Spacing  phantoms 

Visual  Analysis  of  DTS  Spacing  Phantom  Images 

The  350  and  300  pm  spacing  phantoms  could  be  resolved  visually  when  oriented  0°,  22.5°,  and  67.5°  to  the 
scanning  direction.  The  250  and  200  pm  phantoms  could  only  be  partially  resolved,  while  150  and  100  pm 
phantoms  had  minimal  texture  and  no  visually  detectable  banding  pattern.  At  the  90°  orientation,  textural 
information  was  completely  lost  for  all  spacing  phantoms. 

Quantitative  Analysis  of  Microcomputed  Tomography  and  DTS  Spacing  Phantom  Images 
Average  FFT-measured  layer  thickness  agreed  well  with  thicknesses  measured  by  micrometer  during  phantom 
preparation,  supporting  the  use  of  pCT  images  as  an  established  reference  for  further  analyses. 

For  DTS  images,  of  the  four  phantom  angles  relative  to  the  scanning  direction,  orienting  the  length  of  the 
spacing  phantom  in  line  with  the  scanning  direction  yielded  the  greatest  magnitude  of  frequency  (0°  off 
scanning  direction  in  Figure  32a).  Taking  the  0°  angle  as  a  best  resolution  case,  the  FFT  analysis  method  on 
DTS  images  accurately  identified  the  phantom  spacing  for  350  pm  and  300  pm,  but  the  peaks  were  lost  at 
narrower  spacing  (higher  frequencies)  consistent  with  the  Nyquist  criterion  (Figure  32b). 


Figure  32:  (a)  FFT  magnitude  was  highest  when  layers  of  spacing  phantom  were  oriented  perpendicular  to  the 
scanning  direction,  (b)  For  the  0°  angle  scan,  the  FFT  method  was  able  to  identify  spacing  accurately  for  350 
and  300  pm  spacing  phantoms,  but  250  pm  and  below  could  not  be  resolved. 

In  the  MIL,  LFD  and  fractal  analyses,  we  interpret  the  results  with  the  goal  of  finding  parameters  that  are 
sensitive  to  changes  in  spacing,  i.e.,  those  that  monotonically  increase  or  decrease  with  increasing  spacing 
and  those  that  correlate  with  pCT  measurements  of  the  same  kind. 

Most  of  the  MIL  and  LFD  parameters  measured  from  DTS  did  not  follow  a  monotonic  relationship  with  spacing 
for  the  full  range  of  spacing,  although  some  trends  could  be  visible  (Figure  33a,  b).  LFD  parameters  generally 
performed  better  than  MIL  parameters.  Consequently,  higher  correlations  were  found  between  DTS  LFD  and 
pCT  LFD  parameters  than  between  DTS  MIL  and  pCT  MIL  parameters  (R2=0.60  for  MIL.SD  vs  R2=0.64-0.95 
for  LFD. Max,  LFD.EIlipseMax,  LFD. DA,  and  LFD.SD).  Especially  noteworthy  is  the  anisotropy  of  LFD  for  which 
DTS  and  pCT  correlated  95%. 

Fractal  dimension  decreases  with  phantom  spacing  (increasing  image  homogeneity)  until  150  pm,  below  which 
aluminum  is  inseparable  from  plastic  to  resolve  a  scale-dependent  pattern  (Figure  33c).  The  decrease  in 
lacunarity  following  a  peak  at  250  pm  is  possibly  because  texture  heterogeneity  is  no  longer  detected  due  to 
low  resolution  and  the  image  appears  uniform  to  the  lacunarity  analysis  (Figure  33c). 

Thickness  Phantoms 

Visual  Analysis  of  DTS  Thickness  Phantom  Images 

Aluminum  thicknesses  of  120,  115,  90,  75,  and  60  pm  phantoms  could  be  resolved  visually  when  oriented  0°, 
22.5°,  and  67.5°  to  the  scanning  direction,  while  25  pm  phantoms  had  minimal  texture  and  no  visually 
detectable  banding  pattern.  At  the  90°  orientation,  a  banding  pattern  could  be  identified  visually  in  120,  115, 

90,  and  75  pm  phantoms,  but  not  in  60  and  25  pm  phantoms. 
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Figure  33:  The  changes  in  morphometric  parameters  with  varying  phantom  spacing  as  measured  from  DTS 
images.  All  parameters  are  normalized  for  comparison  purposes,  (a)  MIL,  (b)  LFD,  (c)  fractal  and  (d)  FFT 

parameters. 

As  with  the  spacing  phantoms,  thicknesses  were  resolved  best  in  the  0°  scan,  well  in  the  oblique  angles,  and 
worst  in  the  90°  scan.  This  is  demonstrated  as  a  weakened  response  of  FFT  magnitude  at  the  same  frequency 
at  different  specimen  orientations  (Figure  34a).  FFT  magnitude  decreased  together  with  decreasing  layer 
thickness  until  75  pm,  below  which  aluminum  could  not  be  resolved  from  plastic  by  the  analysis  method.  Taken 
together  with  spacing  phantom  results,  FFT  appears  to  be  a  sensitive  tool  for  determining  presence  of  small, 
sufficiently  separated  features  within  a  structure.  The  effect  of  angle  on  resolving  power  emphasizes  the 
importance  of  performing  serial  scans  at  various  angles  to  extract  the  maximum  amount  of  available  structural 
information. 

Although  FFT  analysis  was  able  to  resolve  frequencies  down  to  75  pm,  a  banding  pattern  was  still  visible  in  the 
60  and  25  pm  phantoms,  albeit  faint.  For  this  reason,  MIL,  LFD,  and  fractal  analyses  were  again  performed  on 
regions  cropped  from  pCT  and  DTS  images.  MIL. Max,  MIL. DA,  and  MIL.SD  decreased  with  plate  thickness 
(Figure  35a).  Of  the  aforementioned  MIL  parameters,  MIL. Max  and  MIL.SD  were  moderately  correlated  to  pCT 
MIL  measurements  (R2=0.68  and  0.72,  respectively).  Among  the  LFD  parameters,  LFD. Max,  LFD.SD,  and 
LFD.EIlipseMax  decreased  monotonically  with  thickness  down  to  25  pm.  Of  these  LFD  parameters,  LFD. Max, 
LFD.Av,  LFD.SD,  and  LFD.EIlipseMax  were  highly  correlated  with  pCT  LFD  measurements  (R2=0.97,  0.97, 
0.95,  and  0.98,  respectively).  Fractal  dimension  remained  relatively  constant  while  lacunarity  decreased 
monotonically  with  thickness  down  to  25  pm.  These  findings  suggest  MIL,  LFD  and  fractal  analyses  can  detect 
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changes  in  feature  thickness.  However  these  changes  in  the  measured  parameters  are  of  correlational  nature 
with  respect  to  the  underlying  structural  changes  at  small  resolutions. 


(b) 

Figure  34:  a)  FFT  magnitude  was  highest  when  layers  of  thickness  phantom  were  oriented  perpendicular  to 
the  scanning  direction  and  lower  when  layers  were  parallel,  b)  For  the  0°  angle  scan,  the  FFT  method  was  able 
to  identify  thickness  accurately  for  75-120  pm  thickness  phantoms,  but  60  pm  and  below  could  not  be 

resolved. 


Aluminum  grid  panel  phantom 

All  DTS  images  exhibited  a  varying  degree  of  “halo”  artifacts  adjacent  to  the  aluminum  grid  (Figure  36a,  b). 
The  gray  value  of  the  aluminum  grid  was  shown  as  a  valley  as  it  exhibited  a  higher  attenuation  due  to  its  higher 
density  than  air  (Figure  36c).  The  location  of  this  valley  was  used  to  determine  the  location  of  the  grid  and  the 
amplitude  was  used  to  monitor  signal  variations  in  the  DTS  images. 

All  gray  value  profiles  generated  from  measurements  of  11  slices  /stack  in  four  orientations  (0°,  15°,  22.5°  and 
45)  performed  for  both  the  LI  and  L2  ROIs  are  plotted  in  Figure  37.  Within  each  plot  box  in  Figure  37,  the  11 
gray  value  distributions  within  each  image  slice  stack  overlap  nearly  perfectly  along  the  x-axis  (transverse 
distance  from  the  center  of  image).  These  results  indicate  that  there  is  no  significant  sliding  of  peak  locations 
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between  0  mm  to  11  mm  along  the  z-axis  (i.e.,  depth  of  grid  panel)  for  LI  or  L2  measurements.  Furthermore 
the  results  revealed  that,  while  most  of  the  plot  boxes  show  clear  signals  of  peaks  (gray  value  range:  100-20), 
signals  at  0°  performed  poorly  for  the  L2  orientation.  This  confirms  results  from  spacing  and  thickness 
phantoms  at  the  90°  (parallel)  angle,  indicating  that  DTS  performs  poorly  in  reconstructing  structural  features 
that  are  parallel  to  the  scanning  direction  (also  (Flynn  et  al.,  2007)). 
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Figure  35  continued 


Figure  35:  The  changes  in  morphometric  parameters  with  varying  phantom  thickness  as  measured  from  DTS 
images.  All  parameters  are  normalized  for  comparison  purposes,  (a)  MIL,  (b)  LFD,  (c)  fractal  and  (d)  FFT 

parameters. 


Figure  36:  A  typical  gray  value  profile  distribution  in  the  ROI  shown  in  yellow.  The  DTS  image  has  a 
charateristic  “halo”  (A  and  B)  reflecting  a  sudden  change  in  density  (i.e. ,  air-to-aluminum).  The  halo  has  a 
slightly  higher  gray  value  (i.e.,  whiter).  The  air  is  marked  as  C. 


Table  1  shows  the  mean  and  standard  deviations  of  peak  interval  distances  within  each  ROI.  The  table  also 
confirms  that  there  is  a  minimal  variation  of  mean  interval  distance  between  grids  within  each  ROI  setting  as 
well  as  among  the  eight  orientations.  The  results  indicate  that  there  is  no  measurable  distortion  effect  among 
two  different  measurement  orientations  (LI  vs.  L2)  as  well  as  different  scanning  orientations  (0°,  15°,  22.5°  and 
45°).  This  in  turn  assures  that  the  dimension  data  obtained  in  DTS  is  likely  to  be  reliable  given  that  the  subject 
geometry,  structure  and  density  characteristics  are  comparable  to  the  aluminum  grid  used  here. 


Table  1 :  Mean  and  standard  deviations  of  peak  interval  distances  of  different  ROIs. 


Orientation-ROI 

o 

O 

0°  L2 

15°  LI 

15°  L2 

22.5°  LI 

22.5°  L2 

45°  LI 

45°  L2 

No.  grids  in  ROI  (n) 

20 

17 

19 

17 

17 

17 

18 

16 

Mean  interval  distance 
between  two  peaks  (mm) 

12.5 

12.6 

12.6 

12.6 

12.6 

12.6 

12.6 

12.6 

Standard  deviation 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 
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Figure  37:  Gray  value  profile  plots  of  LI  and  L2  ROIs  from  0°,  15°,  22.5°  and  45°  scanning  orientation  for  the 
grid  phantom.  Note  how  all  peak  locations  are  near  perfectly  aligned  between  different  slices  within  each  ROI. 


Mean  peak  values  of  11  slices  from  each  angular  orientation  (n  =  8)  were  plotted  in  Figure  38.  In  general,  the 
data  indicate  that  the  level  of  signal  strength  (gray  value)  is  affected  by  i)  angular  orientation  of  the  grid  relative 
to  the  scan  direction,  ii)  angle  of  measurements  (LI  vs  L2)  within  each  scanning  orientation  and  iii)  distance 
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from  the  center  of  the  scanning  area. 

While  0°  LI  and  45°  LI  generated  highest  signal  levels,  the  signal  from  their  L2  counterparts  was  much  less;  0° 
L2  generated  nearly  no  signal  and  45°  L2  generated  approximately  40%  less  signal  than  that  of  45°  LI.  The 
differences  between  LI  and  L2  were  reduced  in  15°  and  22.5°  sets.  Overall,  the  22.5°  scan  performed  best  as 
LI  and  L2  signals  were,  as  a  pair,  relatively  uniform  across  the  distance  from  the  center  of  the  image  and  had 
comparable  signal  strength.  The  near-zero  signal  of  0°  L2  confirms  the  previous  observation  that  features  that 
are  parallel  to  the  direction  of  tomographic  acquisition  are  not  well  resolved  in  the  final  DTS  image  (also  (Flynn 
et  al.,  2007)). 

The  signals  varied  as  the  distance  from  the  center  of  the  detector  panel  varied,  for  example  from  the  left  (-ve) 
to  right  (+ve),  the  signals  of  45°  LI,  L2  and  22.5°  LI  increased,  0°  LI  decreased  and  15°  L2,  22.5°  L2  formed  a 
plateau  at  the  center  and  their  signals  decreased  at  the  edges  of  measurement  area.  In  general  scans 
performed  in  oblique  orientations  resulted  in  more  uniform  signals  across  the  ROI. 

Based  on  the  results  of  the  grid  study,  it  is  recommended  that  additional  scans  performed  at  22.5°  orientation 
should  be  explored  for  objects  that  contain  both  vertical  and  horizontal  features  such  as  vertebral  bodies. 
Furthermore  the  finding  that  signal  levels  in  the  22.5°  scanning  orientation  are  relatively  uniform  across  the 
measuring  area  of  30  x  30  cm  (±  150mm  as  shown  in  Figure  38)  supports  this  approach  as  this  ROI  is  large 
enough  to  image  six  to  eight  levels  of  spine  at  once. 


Figure  38:  Mean  peak  gray  value  profiles  of  all  ROI  orientations.  The  plots  are  drawn  upside  down  to  show 
signal  strength  relative  to  the  background  signal  (=  -130  in  blue). 


Results  from  Tasks  3  and  4 


Phase  1  Results 


Relationship  between  pCT  stereology  and  DTS  (LFD,  MIL  and  fractal)  parameters: 


Relationships  were  examined  between  pCT-derived  and  DTS-derived  parameters  of  cancellous  microstructure. 
The  analysis  consisted  of  mixed  linear  regression  models  between  stereological  parameters  of  pCT  (BV/TV, 
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Tb.Th,  Tb.N,  Tb.Sp,  DA)  and  DTS  parameters  (fractal,  MIL  and  LFD).  Mixed  models  were  used  in  lieu  of  simple 
regression  models  to  account  for  pseudoreplication  in  the  data.  Statistical  analysis  was  carried  using  JMP  10 
(Cary,  NC). 

MIL  and  LFD  parameters  calculated  from  the  aggressively  binarized  (Bin.Ag)  images  had  stronger  and  larger 
number  of  significant  relationships  with  pCT  parameters  than  those  from  the  unsharpened  (USM)  and 
conservatively  binarized  (Bin.Csv)  methods  (Table  2).  This  has  resulted  in  no  single  filter  method  (i.e.,  USM, 
Bin.Csv,  Bin.Ag)  producing  a  complete  set  of  DTS-derived  fractal  or  morphometric  parameters  that  correlated 
with  all  of  the  pCT  parameters  tested.  Fractal  dimension  and  lacunarity  slope  parameters  from  fractal  analysis 
also  had  strong  relationships  with  pCT  parameters.  However,  mean  lacunarity  (A)  was  not  related  to  the  pCT 
parameters.  Interestingly,  neither  the  magnitude  nor  the  heterogeneity  of  Tb.Sp  (Tb.Sp. Av  and  Tb.Sp. SD, 
respectively)  and  heterogeneity  of  Tb.N  (Tb.N.SD)  were  related  to  DTS  parameters. 

Table  2:  Coefficients  of  determination  (R2)  from  mixed  linear  models  where  one  pCT  parameter  (columns)  was 
the  outcome  and  one  DTS  parameter  (rows)  was  the  effect  variable.  This  table  was  constructed  using  mixed 
models  instead  of  using  simple  correlations  to  account  for  pseudo-replication.  Direction  of  the  correlation 
coefficient  (f  for  positive  and  [  for  negative  relationships)  and  the  p-value  (in  parenthesis)  are  given  in  each 
cell  along  with  the  R  value.  Only  those  models  with  p  <  0.05  are  shown. 


R* 

(p-value) 

BVTV.Av 

Tb.Th. Av 

Tb.N.Av 

Tb.Sp.Av 

BVTV.SD 

Tb.Th.SD 

Tb.N.SD  Tb.Sp.SD  pCT.DA 

Fractal 

FD 

|0.72 

(<0.002) 

i  0.65 
(0.01) 

A 

Sa 

i  0.89 
(0.002) 

CO  o 

CO  O 

d  P 

— 

i  0.82 
(0.001) 

i  0.77 
(<0.0001 ) 

USM 

(Unsharp 

MIL. Max 

|0.93 

(<0.02) 

|0. 92 
(0.03) 

t  0.55 
(0.03) 

mask) 

MIL. DA 

t0. 58 
(0.01) 

MlL.Av 

|0.93 

(0.006) 

t  0.92 
(0.01) 

MIL.SD 

t  0.57 
(0.01) 

LFD. Max 

t  0.67 
(<0.009) 

LFD.DA 

LFD.Av 

|0.92 

(0.03) 

t  0.66 
(0.009) 

LFD.SD 

t  0.66 
(0.01) 

Bin.Csv 

(Conservative 

MIL. Max 

t  0.63 
(0.03) 

binarization) 

MIL. DA 

t  0.58 
(0.01) 

MlL.Av 

MIL.SD 

t0. 57 
(<0.02) 

LFD.Max 

t  0.88 
(<0.02) 

t  0.76 
(<0.0001 ) 

t  0.72 
(<0.001 ) 

LFD.DA 

LFD.Av 

t  0.88 
(<0.02) 

t  0.74 
(0.0001) 

t  0.77 
(<0.04) 

t  0.71 
(<0.002) 

LFD.SD 

t  0.70 
(<0.002) 

t  0.65 
(<0.02) 

Bin.Ag 

(Aggressive 

MlL.Max 

binarization) 

MIL. DA 

t  0.59 
(<0.03) 

t  0.62 
(0.03) 

t  0.58 
(<0.02) 

MlL.Av 
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Table  2,  continued. 


MIL.SD 

t  0.57 
(<0.02) 

LFD. Max 

t  0.90 
(<0.004) 

t  0.79 
(<0.0001 ) 

t  0.80 
(0.02) 

t  0.66 
(0.01) 

LFD. DA 

LFD.Av 

t  0.88 

t  0.73 

t  0.69 

(0.03) 

(0.0008) 

(0.003) 

LFD.SD 

t  0.89 

t  0.71 

(0.01) 

(0.001) 

Relationship  between  pCT  (LFD,  MIL  and  fractal)  and  DTS  (LFD,  MIL  and  fractal)  parameters: 

To  examine  the  comparability  between  pCT  and  DTS  in  terms  of  LFD-MIL  parameters,  linear  regression 
analyses  were  performed  between  the  MIL  and  LFD  data  extracted  from  DTS  and  pCT  at  the  same  slice  location. 
Because  a  single  DTS  slice  represents  a  1  mm  layer  while  a  slice  from  pCT  represents  0.045  mm,  20  pCT  slices 
were  individually  analyzed  and  then  averaged  to  span  the  thickness  represented  in  a  DTS  slice. 

To  account  for  the  pseudo-replication  effects,  donor  effect  was  treated  as  a  random  effect  in  the  regression 
models  (Table  3).  Only  MIL. Max,  MIL. DA  and  MIL.SD  of  Bin.Csv  and  USM  were  significant  (p=0.006  to 
p=0.0427).  MIL. Max  parameters  had  the  highest  explanatory  power  (R2adj  =  0.58  and  0.591)  followed  by  the 
MIL.SD  parameters  (R2adj=  0.527,  0.535). 

Table  3:  Summary  of  the  relationships  between  pCT-derived  LFD,  MIL  and  fractal  parameters  versus  DTS- 
derived  counterpart.  Mixed  regression  models  were  used  for  the  analysis.  Only  those  with  p  <  0.05  are  shown. 
For  significant  results,  R2adj,  coefficient,  and  p-value  are  reported  below.  Nonsignificant  results  are  indicated  by 
and  analysis  types  not  applicable  to  a  given  variable  are  indicated  by  NA. 


Mixed  univariate  regression 

Y(|iCT)  vs  X  (DTS) 

Bin.Ag 

Bin.Csv 

USM 

No 

processing 

LFD. Max 

- 

- 

- 

NA 

LFD.Av 

- 

- 

- 

NA 

LFD.DA 

- 

- 

- 

NA 

LFD.SD 

- 

- 

- 

NA 

MIL. Max 

- 

0.580 

8.653 

(0.0427) 

0.591 

3.340 

(0.0193) 

NA 

MlL.Av 

- 

- 

- 

NA 

MIL. DA 

- 

0.296 

0.770 

(0.006) 

0.302 

0.541 

(0.0477) 

NA 

MIL.SD 

- 

0.527 

6.911 

(0.0117) 

0.535 

3.584 

(0.0115) 

NA 

FD 

NA 

NA 

NA 

- 

A 

NA 

NA 

NA 

- 

Sa 

NA 

NA 

NA 

- 

While  pCT  and  DTS  are  fundamentally  different  in  terms  of  acquisition  and  reconstruction  methods,  in  this 
study,  we  have  demonstrated  that  by  use  of  the  MIL  method,  there  is  a  good  degree  of  association  between 
the  two  modalities.  LFD  and  fractal  parameters  failed  to  correlate  with  pCT  parameters,  which  may  be  due  to 
the  fact  that  these  two  methods  are  more  sensitive  to  the  grayscale  and  binary  representation  of  bone  and 
non-bone  phases.  (The  grayscale  values  of  nonbone  phases  in  pCT  are  generally  not  useful.)  Other  possible 
reasons  may  be  that  the  fact  that  while  each  DTS  slice  represented  an  integration  of  images  over  a  1  mm  thick 
slice  of  bone,  pCT  data  was  an  average  of  20  individually  analyzed  slices.  This  study  thus  motivates  the  use  of 


39 


MIL  in  DTS  images  as  parameters  representative  of  true  MIL  properties.  However,  the  slopes  of  regressions 
were  greater  than  1  indicating  that  DTS  may  have  precision  but  lacks  accuracy  in  MIL  and  may  require  further 
calibrations.  The  lack  of  an  association  between  DTS  derived  and  pCT  derived  LFD  parameters  but  presence 
of  correlations  between  DTS  derived  LFD  and  pCT  derived  non-LFD  microstructural  parameters  indicate  that 
LFD  measurements  contain  information  on  the  texture  of  cancellous  bone  structure  but  this  information  is  of 
correlational  nature  and  unlikely  a  true  representation  of  the  bone’s  LFD  properties. 

Relationship  between  image  pre-processing  and  correlations  between  DTS  parameters  and  FE-derived 
mechanical  parameters: 

The  ability  of  DTS-derived  parameters  to  predict  bone  stiffness  and  stress  distribution  properties  (Efem,  VMCV 
and  VMExp/aapp),  was  examined  with  and  without  BV/TV.Av  in  the  regression  equations.  BV/TV.Av  was 
introduced  to  the  regression  models  to  examine  predictive  capability  of  DTS  parameters  beyond  that  of  BV/TV.Av 
alone,  in  which  BV/TV.Av  was  treated  as  a  surrogate  measurement  of  bone  mass  (i.e.,  clinical  BMD). 

Statistical  significance  was  set  at  p  <  0.05.  The  multiple  regression  models  were  performed  in  a  stepwise  manner 
by  eliminating  parameters  that  exhibit  multi-collinearity  from  the  initial  predictor  set  based  on  variance  inflation 
factors  (VIF  >  5)  and  then  eliminating  parameters  with  nonsignificant  effects  (p  >  0.05). 

USM  method: 


DTS  images  processed  using  the  USM  method  provided  parameters  that  were  comparable  to  bone  mass  in 
their  ability  to  explain  EFem  and  VMExp/oap  but  more  explanatory  than  BV/TV  alone  for  VMCV  (Table  4).  When 
used  together  with  BV/TV.Av  in  multi  regression  models  (simulating  a  situation  where  bone  mass  information  is 
already  available  from  another  modality),  addition  of  DTS-derived  parameters  contributed  significantly  to 
prediction  of  EFEm  and  VMCV  but  not  VMExp/aapp:  In  addition  to  BVATV,  SA,  LFD. Max,  MlL.Av  and  MIL. DA 
significantly  contributed  to  EFEm  and  FD  and  MIL. DA  to  VMCV.  (Figure  39) 

Table  4:  Summary  of  multiple  regression  models  of  BV/TV-only,  DTS-only  and  DTS+BV/TV  for  the  predictions 
of  Efem,  VMCV,  VMExp/oapp  using  unsharp  mask  method  (USM).  The  table  cells  for  model  terms  include  the 
variable  name,  estimates  of  the  coefficient  and  p-values  associated  with  that  variable.  Only  those  with  p  <  0.05 
are  shown. 


USM  Model 

R'adj 

pCT  Term  DTS  Term  1  DTS  Term  2  DTS  Term  3  DTS  Term  4  Intercept 

Efem  BV/TV 

0.82 

BVTV.Av 

1694  -61.91 

(<0.0001 ) _ (0.05) 


DTS 

0.88 

Sa 

-2994 

(0.0001) 

LFD. Max 

1579 

(0.005) 

MlL.Av 

-59.21 

(0.0009) 

MIL. DA 

247.8 

(0.0004) 

-218.2 

(0.03) 

BV/TV+DTS 

0.93 

BVTV.Av 

SA 

LFD.Max 

MlL.Av 

MIL.DA 

1376 

-1947 

1193 

-30.20 

148.0 

-306.5 

(<0.0001 ) 

(0.0001) 

(0.002) 

(0.005) 

(0.001) 

(0.0001) 

VMCV  BV/TV 

0.57 

BVTV.Av 

-1.3680 

0.9867 

(0.005) 

(0.0001) 

DTS 

0.76 

FD 

LFD.Av 

LFD.SD 

MIL.DA 

1.247 

5.572 

-11.020 

-0.2336 

-2.280 

(0.003) 

(0.03) 

(0.01) 

(0.03) 

(0.04) 

BV/TV+DTS 

0.71 

BVTV.Av 

FD 

LFD.Av 

LFD.SD 

MIL.DA 

-1.277 

0.9519 

- 

- 

-0.2770 

-1.199 

(0.003) 

(0.02) 

- 

- 

(0.01) 

(<0.3) 

VMExp/aapp  BV/TV 

0.96 

BVTV.Av 

-49.34 

14.58 

(0.0001) 

(0.0001) 

DTS 

0.94 

Sa 

MIL. Max 

MIL. DA 

52.23 

0.8621 

-5.748 

14.94 

(0.0001) 

(0.001) 

(0.002) 

(0.0001) 

BV/TV+DTS 

0.96 

BVTV.Av 

Sa 

MIL. Max 

MIL. DA 

-49.34 

- 

- 

- 

14.58 

(0.0001) 

- 

- 

- 

(0.0001) 
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A:  USM  Efem 


B:  USM  VMCV 


Figure  39:  Multiple  regression  results  from  the  unsharp  mask  (USM)  processing  method.  Fit  of  the  multiple 
regression  models  constructed  using  both  BV/TV  and  DTS-based  predictor  variables  and  the  leverage  plots  for 
the  individual  effects  with  the  associated  p-values  for  (A)  EFem  and  (B)  VMCV. 

Bin.Csv  method: 


DTS  images  processed  using  the  Bin.Csv  method  also  provided  parameters  that  were  comparable  to  bone 
mass  in  their  ability  to  explain  EFem  and  VMExp/oap  but  more  explanatory  than  BV/TV  alone  for  VMCV  (Table 
5).  When  used  together  with  BV/TV.Av  in  multiregression  models  (simulating  a  situation  where  bone  mass 
information  is  already  available  from  another  modality),  addition  of  DTS-derived  parameters  contributed 
significantly  to  prediction  of  EFEm  and  VMCV  but  not  VMExp/aapP:  In  addition  to  BV/TV,  and  MIL.DA 
significantly  contributed  to  EFEm  and  VMCV.  (Figure  40) 
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A:  Bin.Csv  EFem 


Predicted  EFEm  BVTV.Av  Sa  MIL.DA 


B:  Bin.Csv  VMCV 


Figure  40:  Multiple  regression  results  from  the  conservatively  binarized  method  (Bin.Csv).  Fit  of  the  multiple 
regression  models  constructed  using  both  BV/TV  and  DTS-based  predictor  variables  and  the  leverage  plots  for 
the  individual  effects  with  the  associated  p-values  for  (A)  EFem  and  (B)  VMCV. 

Table  5:  Summary  of  multiple  regression  models  of  BV/TV-only,  DTS-only  and  DTS+BV/TV  for  predictions  of 
EFem,  VMCV,  VMExp/aapp  using  conservatively  binarized  method  (Bin.Csv).  The  table  cells  for  model  terms 
include  the  variable  name,  estimates  of  the  coefficient  and  p-values  associated  with  that  variable.  Only  those 
with  p  <  0.05  are  shown. 


Bin.Csv  Model 

R'adj 

pCT  Term  DTS  Term  1  DTS  Term  2  Intercept 

Efem  BV/TV 

0.82 

BVTV.Av 

1694  -61.91 

(<0.0001 ) _ (0.05) 


DTS 

0.84 

Sa 

-3294 
(<0.0001 ) 

MIL.DA 

220.9 

(<0.004) 

-371.7 

(<0.002) 

BV/TV+DTS 

0.92 

BVTV.Av 

1534 

(<0.0001 ) 

Sa 

-2028 
(<0.0001 ) 

MIL.DA 

175.0 

(<0.0003) 

-433.1 
(<0.0001 ) 

VMCV 

BV/TV 

0.57 

BVTV.Av 

-1.368 

(0.005) 

0.9867 
(<0.0001 ) 

DTS 

0.76 

Sa 

MIL.DA 

3.572 

-0.4306 

1.640 

(<0.0002) 

(<0.0003) 

(<0.0001 ) 

BV/TV+DTS 

0.75 

BVTV.Av 

Sa 

MIL.DA 

-1.015 

2.815 

-0.4041 

1.692 

(<0.02) 

(0.003) 

(<0.0004) 

(<0.0001 ) 

VMExp/CJapp 

BV/TV 

0.96 

BVTV.Av 

-49.34 

14.58 

(<0.0001 ) 

(<0.0001 ) 

DTS 

0.91 

Sa 

46.56 

11.16 

(0.001) 

(<0.0001 ) 

BV/TV+DTS 

0.96 

BVTV.Av 

Sa 

-49.34 

- 

14.58 

(<0.0001 ) 

- 

(<0.0001 ) 
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Bin.Aq  method: 


DTS  images  processed  using  the  Bin.Ag  method,  similar  to  other  methods  above,  provided  parameters  that 
were  comparable  to  bone  mass  in  their  ability  to  explain  EFem  and  VMExp/oap  but  more  explanatory  than  BV/TV 
alone  for  VMCV  (Table  6).  When  used  together  with  BV/TV. Av  in  multiregression  models  (simulating  a 
situation  where  bone  mass  information  is  already  available  from  another  modality),  addition  of  DTS-derived 
parameters  contributed  significantly  to  prediction  of  EFem  and  VMCV  but  not  VMExp/oapp:  In  addition  to  BV/TV, 
and  LFD.SD  significantly  contributed  to  EFEm  and,  MIL. DA  to  VMCV.  In  this  case,  there  seems  to  be  a  large 
confounding  effect  of  BV/TV  on  most  significant  DTS  parameters  such  that  inclusion  of  BV/TV  is  significantly 
worse  than  DTS  alone.  (Figure  41) 


Table  6:  Summary  of  multiple  regression  models  of  BV/TV-only,  DTS-only  and  DTS+BV/TV  for  predictions  of 
EFem,  VMCV,  VMExp/oapp  using  aggressively  binarized  method  (Bin.Ag).  The  table  cells  for  model  terms 
include  the  variable  name,  estimates  of  the  coefficient  and  p-values  associated  with  that  variable.  Only  those 
with  p  <  0.05  are  shown. 


Bin.Ag  Model 

R*adj 

|iCT  Term  DTS  Term  1  DTS  Term  2  DTS  Term  3  Intercept 

Efem  BV/TV 

0.82 

BVTV.Av 

1694  -61.91 

(<0.0001 ) _ (0.05) 


DTS 

0.84 

Sa 

LFD.SD 

-2443 

5853 

-116.9 

(0.0003) 

(0.009) 

(0.03) 

BV/TV+DTS 

0.89 

BVTV.Av  SA 

LFD.SD 

1580  -1470 

3194 

-222.4 

(<0.0001 ) 

(0.002) 

(0.009) 

(<0.0001 ) 

VMCV 

BV/TV 

0.57 

BVTV.Av 

-1.368 

0.9867 

0.005 

(<0.0001 ) 

DTS 

0.76 

FD 

MIL.DA 

MIL.SD 

0.9037 

-0.7169 

0.5414 

-0.8906 

(<0.04) 

(0.006) 

(<0.03) 

0.5 

BV/TV+DTS 

0.68 

BVTV.Av 

FD 

MIL.DA 

MIL.SD 

-1.404 

- 

-0.2951 

- 

1.396 

(0.001) 

- 

(0.001) 

- 

(<0.0001 ) 

VMExp/CJapp 

BV/TV 

0.96 

BVTV.Av 

-49.34 

14.58 

(<0.0001 ) 

(<0.0001 ) 

DTS 

0.93 

Sa 

LFD.DA 

LFD.SD 

30.88 

1.115 

-144.8 

10.06 

(<0.03) 

(0.01) 

(0.01) 

(<0.0001 ) 

BV/TV+DTS 

0.96 

BVTV.Av 

Sa 

LFD.DA 

LFD.SD 

-49.34 

- 

- 

- 

14.58 

(<0.0001 ) 

- 

- 

- 

(<0.0001 ) 

Even  though  it  appears  that  DTS  parameters  calculated  using  the  aggressive  binarization  (Bin.Ag)  method  had 
stronger  and  larger  number  of  significant  relationships  with  pCT  parameters,  DTS  parameters  calculated  using 
the  USM  processing  resulted  in  slightly  better  models  of  mechanical  parameters  than  the  other  two  processing 
methods.  This  may  be  because  the  two  binarization  methods  produced  parameters  that  are  largely  correlated 
to  bone  mass  whereas  the  USM  method  produced  parameters  that  contain  heterogeneity  information 
independent  from  bone  mass  (Table  2).  The  USM  method  is  also  a  “grayscale”  method  which  may  correspond 
to  the  gray-level-based  heterogeneous  FE  models  better  than  the  binarized  images  do. 


43 


A:  Bin.Ag  Efem 


B:  Bin.Ag  VMCV 


Figure  41:  Multiple  regression  results  from  the  aggressively  binarized  method  (Bin.Ag).  Fit  of  the  multiple 
regression  models  constructed  using  both  BV/TV  and  DTS-based  predictor  variables  and  the  leverage  plots  for 
the  individual  effects  with  the  associated  p-values  for  (A)  EFEM  and  (B)  VMCV. 


Phase  2  Results 


Relationship  between  pCT  stereology  and  stack-based  DTS  (LFD,  MIL  and  fractal)  parameters  with 
varying  scan  orientation: 

Mean  (and  standard  deviation  where  appropriate)  pCT  stereology  parameters  of  bone  volume  fraction  (BV/TV), 
trabecular  thickness  (Tb.Th),  number  (Tb.N),  separation  (Tb.Sp),  and  degree  of  anisotropy  (DA),  connectivity 
(Conn.Dn),  fractal  dimension  (FD),  structure  model  index  (SMI)  were  correlated  to  stack-based  DTS 
parameters  (i.e.,  mean  and  standard  deviations  of  MIL,  LFD,  and  fractal  parameters)  calculated  from  six 
scanning  orientations  (AP  and  LM  at  angular  (22.5°),  parallel  (0°)  and  perpendicular  (90°))  using  mixed 
regression  models  (Table  7). 

Table  7:  Summary  of  mixed  regression  models  examining  correlations  between  pCT  stereology  parameters 
and  DTS  parameters.  The  DTS  measurements  used  in  these  analyses  were  from  the  cuboidal  VOIs  as  this 
type  of  VOI  is  most  comparable  to  that  used  in  pCT  analyses.  Only  those  with  p  <  0.05  are  shown.  Models 
which  are  both  significant  in  AP  and  LM  are  highlighted.  AP  =  anterior  posterior  (coronal)  view,  LM  =  Lateral- 
medial  (sagittal)  view.  Scan  directions  ang  =  angular  (22.5°),  par  =  parallel  (0°),  perp  =  perpendicular  (90°)  are 
with  respect  to  the  axis  of  the  spine. 


Y  =  pCT 

DTS 

Mean  (averaged  within  stack) 

SD  (intra-heterogeneity  within 
stack) 

Stereology 

Parameter 

View 

Orientation 

Parameter 

R  adj 

Coefficient 

sign 

p-value 

R  adj 

Coefficient 

sign 

p-value 

BVTV.Av 

AP 

ang 

MlL.Av 

0.881 

+ 

0.0067 

MIL. Max 

0.883 

+ 

0.0048 

MIL.SD 

0.867 

+ 

0.0422 

LM 

perp 

LFD.Av 

0.870 

+ 

0.0322 

LFD. Max 

0.874 

+ 

0.0215 

SA 

0.867 

+ 

0.0359 

44 


Table  7,  continued. 


BVTV.SD 

AP 

ang 

MlL.Av 

0.757 

+ 

0.0464 

MIL. Max 

0.758 

+ 

0.0436 

LM 

par 

FD 

0.760 

+ 

0.0363 

LFD.Max 

0.771 

+ 

0.0195 

Tb.Th.Av 

AP 

ang 

MlL.Av 

0.636 

+ 

0.0237 

MIL. DA 

0.621 

+ 

0.0334 

MIL. Max 

0.647 

+ 

0.0155 

MIL.SD 

0.661 

+ 

0.0050 

perp 

LFD.Av 

0.634 

+ 

0.0321 

LFD.Max 

0.659 

+ 

0.0096 

LFD.SD 

0.618 

+ 

0.0462 

MIL. DA 

0.687 

+ 

0.0016 

0.558 

+ 

0.0291 

MIL.SD 

0.702 

+ 

0.0026 

LM 

ang 

FD 

0.581 

+ 

0.0480 

MIL. DA 

0.639 

+ 

0.0046 

par 

MIL. DA 

0.624 

- 

0.0011 

MIL.SD 

0.566 

- 

0.0440 

perp 

LFD.Av 

0.691 

+ 

0.0020 

LFD.Max 

0.705 

+ 

0.0012 

LFD.SD 

0.656 

+ 

0.0118 

MIL. DA 

0.705 

+ 

0.0006 

MIL.SD 

0.704 

+ 

0.0022 

Tb.Th.SD 

AP 

ang 

SA 

0.629 

- 

0.0229 

par 

LFD.DA 

0.623 

+ 

0.0360 

SA 

0.636 

+ 

0.0201 

perp 

LFD.Av 

0.635 

+ 

0.0286 

LFD.Max 

0.680 

+ 

0.0033 

LFD.SD 

0.648 

+ 

0.0117 

MIL. DA 

0.743 

+ 

0.0000 

0.609 

+ 

0.0069 

MIL. Max 

0.656 

+ 

0.0184 

MIL.SD 

0.764 

+ 

0.0000 

0.667 

+ 

0.0029 

LM 

ang 

MIL. DA 

0.600 

+ 

0.0202 

0.603 

+ 

0.0362 

SA 

0.626 

- 

0.0375 

par 

LFD.DA 

0.644 

- 

0.0192 

MIL. DA 

0.627 

- 

0.0277 

SA 

0.613 

- 

0.0394 

perp 

LFD.Av 

0.656 

+ 

0.0238 

LFD.Max 

0.643 

+ 

0.0348 

MIL. DA 

0.772 

+ 

0.0000 

MIL.SD 

0.737 

+ 

0.0003 

Tb.N.Av 

AP 

ang 

Sa 

0.923 

+ 

0.0210 

0.929 

+ 

0.0096 

par 

Sa 

0.927 

+ 

0.0100 

0.920 

+ 

0.0499 

LM 

ang 

Sa 

0.931 

+ 

0.0044 

par 

Sa 

0.930 

+ 

0.0044 

Tb.N.SD 

AP 

ang 

LFD.DA 

0.735 

+ 

0.0309 

LM 

ang 

LFD.Max 

0.762 

+ 

0.0173 

LFD.SD 

0.768 

+ 

0.0155 

par 

LFD.Max 

0.754 

+ 

0.0206 

LFD.SD 

0.753 

+ 

0.0213 

Tb.Sp.Av 

AP 

ang 

LFD.Av 

0.924 

+ 

0.0173 

Sa 

0.927 

+ 

0.0085 

par 

Sa 

0.920 

- 

0.0207 

LM 

ang 

Sa 

0.917 

- 

0.0439 

par 

MlL.Av 

0.917 

+ 

0.0160 

MIL. Max 

0.916 

+ 

0.0200 

Tb.Sp.SD 

LM 

ang 

MlL.Av 

0.870 

+ 

0.0399 

Sa 

0.878 

+ 

0.0106 

par 

MlL.Av 

0.869 

+ 

0.0359 

0.869 

+ 

0.0270 

MIL. Max 

0.865 

+ 

0.0388 

perp 

MIL. Max 

0.875 

+ 

0.0186 

MIL.SD 

0.884 

+ 

0.0065 

Conn.Dn. 

AP 

ang 

FD 

0.762 

- 

0.0203 

Sa 

0.769 

+ 

0.0065 

0.756 

+ 

0.0457 

par 

FD 

0.739 

+ 

0.0439 

LFD.DA 

0.742 

+ 

0.0485 

LFD.Max 

0.749 

+ 

0.0232 

MIL. Max 

0.733 

+ 

0.0443 

MIL.SD 

0.746 

+ 

0.0164 

Sa 

0.792 

+ 

0.0021 

perp 

MIL. DA 

0.770 

- 

0.0049 

MIL. Max 

0.741 

- 

0.0380 

MIL.SD 

0.773 

- 

0.0036 

LM 

ang 

Sa 

0.809 

+ 

0.0006 

0.749 

+ 

0.0319 

Table  7,  continued. 


par 

MIL. DA 

0.771 

+ 

0.0104 

SA 

0.791 

+ 

0.0011 

perp 

LFD.DA 

0.753 

- 

0.0430 

LFD.SD 

0.754 

+ 

0.0458 

MIL. DA 

0.778 

- 

0.0053 

MIL. Max 

0.769 

- 

0.0107 

MIL.SD 

0.801 

- 

0.0009 

SA 

0.738 

+ 

0.0477 

FD.Av 

AP 

ang 

MIL.Av 

0.893 

+ 

0.0138 

MIL. Max 

0.893 

+ 

0.0134 

LM 

ang 

MIL. DA 

0.887 

+ 

0.0289 

perp 

LFD.Av 

0.887 

+ 

0.0262 

LFD.Max 

0.890 

+ 

0.0194 

SA 

0.885 

+ 

0.0413 

FD.SD 

AP 

ang 

A 

0.764 

+ 

0.0048 

0.748 

+ 

0.0105 

par 

A 

0.743 

+ 

0.0246 

0.725 

+ 

0.0441 

MIL. DA 

0.745 

+ 

0.0121 

MIL. Max 

0.754 

+ 

0.0424 

MIL.SD 

0.763 

+ 

0.0090 

perp 

A 

0.711 

+ 

0.0192 

LFD.DA 

0.754 

- 

0.0147 

0.768 

+ 

0.0093 

LFD.Max 

0.717 

- 

0.0228 

LFD.SD 

0.726 

- 

0.0233 

SA 

0.738 

+ 

0.0144 

LM 

ang 

FD 

0.766 

+ 

0.0011 

A 

0.751 

+ 

0.0182 

0.757 

+ 

0.0047 

LFD.Av 

0.754 

+ 

0.0273 

par 

FD 

0.727 

+ 

0.0170 

A 

0.773 

+ 

0.0040 

0.761 

+ 

0.0027 

MIL. DA 

0.846 

+ 

0.0000 

0.787 

+ 

0.0005 

MIL.SD 

0.771 

+ 

0.0056 

0.761 

+ 

0.0040 

Sa 

0.713 

+ 

0.0390 

perp 

FD 

0.715 

+ 

0.0380 

LFD.Av 

0.699 

- 

0.0130 

LFD.Max 

0.708 

- 

0.0061 

LFD.SD 

0.717 

- 

0.0352 

MIL.Av 

0.749 

- 

0.0106 

MIL. DA 

0.735 

- 

0.0430 

MIL. Max 

0.765 

- 

0.0031 

MIL.SD 

0.773 

- 

0.0015 

Sa 

0.739 

+ 

0.0105 

mean  MIL 

AP 

ang 

MIL.Av 

0.648 

+ 

0.0002 

MIL. Max 

0.647 

+ 

0.0002 

MIL.SD 

0.540 

+ 

0.0066 

SD  MIL 

AP 

perp 

Sa 

0.402 

+ 

0.0302 

LM 

ang 

LFD.Av 

0.411 

+ 

0.0341 

LFD.Max 

0.404 

+ 

0.0361 

LFD.SD 

0.400 

+ 

0.0269 

par 

LFD.Av 

0.520 

+ 

0.0041 

0.346 

+ 

0.0234 

LFD.DA 

0.346 

+ 

0.0445 

0.394 

+ 

0.0178 

LFD.Max 

0.479 

+ 

0.0092 

LFD.SD 

0.484 

+ 

0.0036 

perp 

Sa 

0.398 

+ 

0.0283 

SMI 

LM 

perp 

LFD.Av 

0.943 

_ - _ 

0.0151 

LFD.DA 

0.940 

- 

0.0394 

LFD.Max 

0.945 

- 

0.0094 

LFD.SD 

0.945 

- 

0.0096 

A  number  of  significant  correlations  were  found  between  pCT  and  DTS  variables.  In  general,  perpendicular 
scans  in  the  LM  view  provided  significant  correlates  for  a  greater  number  of  pCT  parameter  than  other  scans. 
As  suggested  by  the  phantom  studies  of  Task  2,  it  is  possible  that  the  perpendicular  scans  are  better  for 
examining  the  vertically  aligned  trabeculae  while  parallel  scans  are  better  for  examining  the  horizontally 
aligned  trabeculae.  While  the  perpendicular  scans  in  an  LM  view  appear  to  be  the  preferable  single  scan 
configuration,  it  is  likely  that  scans  in  multiple  orientations  will  provide  complementary  information.  This  is  most 
notable  for  the  cases  of  Tb.N.Av,  Tb.Sp.Av  and  Tb.Sp.SD  where  parallel  scans  provided  significant  correlates 
but  perpendicular  scans  did  not.  Overall,  the  results  support  the  idea  that  DTS  can  provide  quantitative 
information  on  the  cancellous  bone  microstructure  and  possibly  can  contribute  towards  the  prediction  of  bone 
strength. 
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Relationship  between  vertebral  fracture  strength  and  DTS  (LFD,  MIL  and  fractal)  parameters  with 
varying  scan  orientation  and  VOI  method: 


We  found  several  significant  correlations  between  vertebral  fracture  strength  and  DTS-derived  cancellous  bone 
parameters  (Table  8).  The  adaptive  VOI  method,  which  allows  for  examination  of  a  larger  volume  of  cancellous 
bone,  produced  a  higher  number  of  significant  models  with  1  to  3%  higher  explained  variability  compared  to 
the  other  VOI  methods.  Overall,  analysis  of  the  image  stack  provided  more  explanatory  information  than  the 
analysis  of  a  single  slice. 

Table  8:  Summary  of  significant  (p<0.05)  relationships  between  peak  fracture  force  and  (A)  mean(DTS)  and  (B) 
SD(DTS)  parameters  as  assessed  by  mixed  linear  regression.  The  results  are  grouped  by  slice/stack,  ROI-VOI 
methods  and  scanning  orientations.  For  ease  of  interpretation,  the  pCT  parameters  that  are  associated  with 
the  DTS  parameter  under  question  are  presented  in  the  “Correlates  to  pCT”  column.  AP  =  anterior  posterior 
coronal  plane,  LM  =  Lateral-medial  sagittal  plane,  ang  =  angular  (22.5°),  par  =  parallel  (0°),  perp  = 
perpendicular  (90°). _ _ _ _ 


(A)  Mean  (DTS) 

Slice 

Stack 

Largest  ROI 

Cuboidal  VOI 

Adaptive  VOI 

View 

Orientation 

DTS 

Parameter 

Correlates 
to  pCT 

R  adj 

Coeff. 

P- 

value 

R  adj 

Coeff. 

P- 

value 

R  adj 

Coeff. 

P- 

value 

AP 

ang 

Sa 

Tb.Th.SD, 

Tb.N.Av, 

Conn.Dn 

0.415 

30530 

0.050 

perp 

LFD. Max 

FD.SD 

0.456 

-66932 

0.044 

LFD.SD 

FD.SD 

0.453 

293629 

0.032 

MIL. DA 

Tb.Th.Av, 

Conn.Dn. 

0.461 

-2423 

0.048 

MIL. Max 

Conn.Dn. 

0.467 

-606 

0.030 

0.478 

-638 

0.019 

MIL.SD 

Conn.Dn. 

0.504 

-3191 

0.014 

0.514 

-3471 

0.009 

LM 

ang 

Sa 

Tb.Th.SD, 

Tb.N.Av, 

Conn.Dn 

0.437 

31260 

0.024 

par 

Sa 

Tb.N.Av 

0.389 

30278 

0.040 

perp 

LFD.Av 

Tb.Th.Av 

0.509 

129171 

0.009 

0.459 

137404 

0.041 

0.489 

153525 

0.020 

LFD. Max 

FD.SD 

0.475 

-77805 

0.031 

0.490 

-85184 

0.026 

LFD.SD 

FD.SD 

0.455 

275525 

0.024 

0.448 

283332 

0.028 

MIL. DA 

Tb.Th.Av, 

Conn.Dn. 

0.479 

-3485 

0.044 

0.484 

-3634 

0.041 

MIL. Max 

Conn.Dn. 

0.427 

-429 

0.025 

0.458 

-496 

0.011 

MIL.SD 

Conn.Dn. 

0.490 

-2949 

0.007 

0.511 

-3050 

0.005 

(B)  SD  (DTS) 

Stack 

Cuboidal  VOI 

Adaptive  VOI 

View 

Orientation 

DTS 

Parameter 

Correlates 
to  pCT 

R  adj 

Coeff. 

P- 

value 

R  adj 

Coeff. 

P- 

value 

AP 

perp 

LFD. Max 

FD.SD 

0.444 

151319 

0.043 

LFD.SD 

FD.SD 

0.477 

633858 

0.016 

MIL. Max 

Conn.Dn. 

0.458 

-990 

0.029 

LM 

perp 

MIL. Max 

Conn.Dn. 

0.321 

-895.212 

0.037 

0.316 

-1561 

0.007 

MIL.SD 

Conn.Dn. 

0.276 

6436.189 

0.032 

0.330 

-7256 

0.031 
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Multiple  regression  models  constructed  using  DTS  parameters  from  the  adaptive  VOI  approach  indicated  that 
using  a  combination  of  variables  measured  from  a  single  DTS  configuration,  the  prediction  strength  for  vertebral 
fracture  force  can  be  substantially  increased  (up  to  86%)  (Table  9). 

Models  from  the  perpendicular  scanning  orientations  had  greater  explanatory  capability  for  vertebral  strength 
than  did  the  parallel  or  angular  scans.  Parameters  of  vertebral  size  (width  and  height)  were  significant  and 
addition  of  a  vertebral  size  parameter  in  the  model  increased  explanatory  capability  of  the  model  for  six 
scanning  configurations.  The  model  with  the  highest  explanatory  capability  was  for  the  AP/perpendicular 
configuration.  The  addition  of  vertebral  size  into  the  models  changed  the  composition  of  predictors  in  the  model 
for  most  cases.  However,  three  out  of  the  four  predictor  parameters  persisted  up  on  addition  of  vertebral  size  in 
to  the  AP/perpendicular  model,  indicating  some  level  of  robustness  for  these  predictors. 

The  study  thus  demonstrates  mixed  multiple  regression  models  of  DTS  parameters  can  predict  fracture 
strength  of  vertebral  bone  where  the  best  scanning  orientation  is  anterior-posterior  with  DTS  acquisition  motion 
perpendicular  to  the  torso.  We  have  observed  that  vertebral  size  itself  is  a  significant  contributor  to  fracture 
strength  and  further  improves  the  predictability  by  5  to  28%. 


Table  9:  Summary  of  the  mixed  multiple  regression  models  for  peak  vertebral  fracture  force  with  and  without 
height  and  width  data  measured  from  DTS  images.  Each  cell  shows  R2adj  followed  by  the  predictor(s)  in  the 
model  (p-value).  Only  those  with  p  <  0.05  are  shown. _ _ 


View 

Orientation 

DTS  parameters  only 

DTS  parameters  +  size 

AP 

Angular 

0.653 

0.755 

(Anterior-posterior 
coronal  plane) 

(23g) 

Av(F.Sa)  (0.0168) 

AP  Height  (0.0001) 

SD(LFD.SD)  (0.0018) 

SD(LFD.SD)  (0.0002) 

SD(F.Sa)  (0.0260) 

Av(MIL.SD)  (0.0405) 

Av(MIL.Max)  (0.0357) 

Parallel 

0.503 

0.712 

(0s) 

SD(MIL.DA)  (0.0032) 

AP  Height  (0.0001) 

Av(MIL.DA)  (0.0431) 

SD(MIL.DA)  (0.0102) 

Perpendicular 

(90g) 

0.814* 

0.860* 

Av(MIL.SD)  (<0.0001 ) 

AP  Height  (0.0006) 

SD(MIL.DA)  (<0.0001 ) 

Av(MIL.SD)  (0.0001) 

Av(MIL.Av)  (0.0028) 

SD(MIL.DA)  (0.0001) 

SD(LFD.SD)  (0.0090) 

Av(MIL.Av)  (0.0015) 

SD(MIL.Max)  (0.0346) 

LM 

Angular 

0.437 

0.669 

(Lateral-medial 
sagittal  plane) 

(23g) 

Av(F.Sa)(0.0243) 

LM  Height  (0.0001) 

Parallel 

0.389 

0.669 

(0s) 

Av(F.Sa)  (0.0398) 

LM  Height  (<0.0001 ) 

Perpendicular 

(90g) 

0.511 

0.713 

Av(MIL.SD)  (0.0045) 

LM  Height  (0.0001) 

SD(LFD.Av)  (0.0141) 

KEY  RESEARCH  ACCOMPLISHMENTS 

We  produced  strong  evidence  that  the  “concept”  being  explored,  i.e.,  that  quantitative  analysis  of  bone  quality 
using  DTS  is  feasible. 

We  characterized  the  DTS  imaging  system  and  identified  its  ability  to  resolve  features  under  different 
operational  parameters. 

We  determined  the  potential  of  the  DTS  system  for  quantification  of  density. 

We  established  methods  to  be  used  for  quantitative  analysis  from  DTS  images  of  vertebral  bone. 
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To  the  best  of  our  knowledge,  this  is  the  first  demonstration  that  correlates  of  conventionally  used  parameters 
of  cancellous  bone  microstructure  can  be  derived  from  quantitative  analysis  of  bone  texture  from  DTS  images. 


We  demonstrated  that  parameters  of  cancellous  bone  microstructure  derived  from  DTS  increases  accuracy  in 
prediction  of  cancellous  bone  stiffness. 

We  demonstrated  that  parameters  of  cancellous  bone  microstructure  and  vertebral  geometry  derived  from  DTS 
increases  accuracy  in  prediction  of  vertebral  strength. 

We  have  identified  that  AP-perpendicular  scanning  orientation  may  provide  optimum  DTS  parameters  for 
prediction  of  bone  fracture  strength. 

We  designed  and  fabricated  an  innovative,  cost-effective  apparatus  that  can  apply  multiaxial  “wedge”  loads 
using  the  capability  of  a  conventional  uniaxial  materials  testing  system. 
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Manuscript  in  preparation  to  be  submitted  to  Spine.  Final  draft  is  under  review  by  the  PI. 
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Bone:  The  Effect  of  Positioning  and  Scan  Orientation  on  Prediction  of  Cancellous  Bone  Stiffness." 

Conference  abstract  submitted  to  the  60th  Annual  Meeting,  Orthopaedic  Research  Society,  March  15-18, 
2014,  New  Orleans,  Louisiana. 

Kim,  W.,  Oravec,  Divine,  G.  W.,  Maatman,  T,  Flynn,  M.  J.  and  Yeni,  Y.  N.,  "Digital  Tomosynthesis  for 
Prediction  of  Human  Whole  Vertebral  Stiffness" 

Conference  abstract  submitted  to  the  6Cfh  Annual  Meeting,  Orthopaedic  Research  Society,  March  15-18, 
2014,  New  Orleans,  Louisiana. 

Presentations 

Structural  and  Microstructural  Considerations  in  Vertebral  Bone  Fragility,  Wayne  State  University,  Department 
of  Biomedical  Engineering,  February  23,  2012. 

Invited  seminar  as  a  part  of  a  graduate  class:  BME  seminar  class  (BME  8070). 

Deciphering  trabecular  microstructure  using  digital  tomosynthesis  (DTS),  Bone  and  Joint  Center  Seminar 
Series,  Henry  Ford  Hospital,  August  3,  2012. 

Institutional  level  seminar. 

Assessing  fracture  risk  in  osteoporotic  spine  using  DTS  -  scanning  orientation  effects  explored,  Bone  and  Joint 
Center  Seminar  Series,  Henry  Ford  Hospital,  June  7,  2013. 

Institutional  level  seminar. 
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Funding  applied  for  based  on  work  supported  by  this  award 

April,  2012:  “Novel  Digital  Tomosynthesis  and  Noncoding  RNA  Approaches  to  Predict  Progression  of 
Osteoarthritis” 

FY12  POD  CDMRP  Peer  Reviewed  Medical  Research  Program  Investigator  Initiated  Award 
Role  on  grant:  Principal  Investigator  (MPI:  Gibson) 

The  proposed  project  was  based  on  the  methods  developed  for  analysis  of  bone  quality  from  DTS  images 
and  the  encouraging  results  obtained  from  vertebral  bone.  However,  we  were  not  invited  to  submit  a  full 
application. 

Tomosynthesis  for  assessment  of  hip  bone  quality  and  fracture  strength 
National  Institutes  of  Health  (R21) 

Requested  Award  period:  1/1/2014  -  12/31/2015 
Role  on  grant:  Principal  Investigator 

This  is  an  exploratory  project  to  examine  the  feasibility  of  DTS  for  assessment  of  hip  bone  quality  and 
fracture  risk.  The  methods  developed  for  analysis  of  bone  quality  from  DTS  images  and  the  encouraging 
results  obtained  from  vertebral  bone  during  the  current  award  support  feasibility  of  the  proposed  work. 
However,  as  we  learned  during  this  award,  scan  configuration  and  analysis  methods  can  be  critical  for  a 
successful  application  to  the  hip  and  these  may  not  be  the  same  as  for  the  spine.  Therefore,  there  is  merit 
in  conducting  research  to  explore  feasibility  and  optimize  methodologies  for  this  problem.  This  proposal  is 
currently  pending  review. 

Tomosynthesis  of  the  human  spine  in  vivo  for  assessment  of  vertebral  bone  quality 
National  Institutes  of  Health  (R01) 

Requested  Award  period:  4/1/2014  -  3/31/2019 
Role  on  grant:  Principal  Investigator 

Having  established  feasibility  and  developed  methods  to  measure  vertebral  bone,  we  prepared  this 
proposal  to  examine  the  in  vivo  feasibility  of  the  approach.  We  will  examine  subjects  with  a  prevalent 
vertebral  fracture  and  subjects  with  primary  hyperparathyroidism  (pHPT)  as  test  cases  representing 
increased  fracture  risk  but  with  different  bone  microstructural/structural  factors  involved.  These  groups  will 
be  compared  to  normal  controls  (without  fracture  or  pHPT)  and  the  ability  of  DTS  to  discriminate  among 
these  groups  will  be  examined  and  compared  to  BMD  from  DEXA,  serum  and  urine  biomarkers  of  bone 
turnover,  bone  remodeling  dynamics  as  determined  from  iliac  biopsies  and  high  resolution  microstructural 
analysis  of  iliac  bone.  This  proposal  will  be  submitted  by  November  12,  2013. 

Employment  or  research  opportunities  applied  for  and/or  received  based  on  experience/training 
supported  by  this  award 

This  project  provided  training  opportunities  for  the  following  individuals  whether  or  not  they  were  compensated 
by  the  project  funds: 

1.  Woong  Kim,  PhD:  University  of  Auckland,  New  Zealand,  Post-doc  fellow,  2012  -  present. 

Scanning  and  analysis  methods  for  DTS  imaging  of  vertebral  bone. 

2.  Richard  Banglmaier,  PhD:  Wayne  State  University,  Senior  Research  Engineer,  2011  -  April  2013. 
Nonuniform  mechanical  testing  of  vertebral  bone. 

3.  Daniel  Oravec,  MSc:  Tampere  University  of  Technology,  Finland,  Research  Engineer,  2011  -  present 
All  aspects  of  the  project. 

4.  Ryan  Bylsma,  Student,  Wayne  State  University,  School  of  Medicine,  October  201 1  -  January  2012. 

Micro-CT  scanning  of  bone 

5.  Mary  Nixon,  Student,  Wayne  State  University,  School  of  Medicine,  May  2012  -  July  2012. 

Image  processing  and  analysis. 

6.  Justin  Schupbach,  Student,  Wayne  State  University,  School  of  Medicine,  June  2012  -  July  2012. 

Image  processing  and  analysis. 
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7.  Kaitlin  McLoughlin,  Student,  Wayne  State  University,  School  of  Medicine,  August  2012  -  September  2012. 
Stereologic  analysis  of  cancellous  bone. 

8.  Kalyan  Sreeram,  Student,  Wayne  State  University,  School  of  Medicine,  August  2012  -  April  2013. 

Image  processing  and  analysis. 

9.  Thomas  K  Maatman,  Student,  Wayne  State  University,  School  of  Medicine,  August  2012  -  May  2013. 
Micro-CT  scanning  of  bone,  bone  phantoms  and  image  analysis.  Finite  element  modeling. 

10.  Nicole  Ramo,  Student,  Kettering  University,  Biomedical  Engr.,  Winter  2011,  Spring  2012. 

Micro-CT  image  analysis.  Specimen  preparation. 

11.  Jason  Bagnall,  Student,  Oakland  Community  College,  Engineering,  March  -  August  2012. 

Micro-CT  image  processing.  Stereologic  analysis  of  cancellous  bone. 

12.  Vincent  Giacopelli,  Student,  Harrison  High  School,  Farmington,  Ml,  June  18  -  August  10,  2012. 

DTS  image  processing,  LFD  analysis  of  cancellous  bone  images. 

13.  Toufic  Jildeh,  Student,  Wayne  State  University,  School  of  Medicine,  May  2013  -  September  2013. 

DTS  and  radiographic  image  processing  and  analysis. 

14.  Joy  Zhang,  Student,  Detroit  Country  Day  High  School,  Beverly  Hills,  Ml,  June  10  -  August  9,  2013. 
Micro-CT  image  processing.  Stereologic  analysis  of  cancellous  bone. 

15.  Angela  Xiao,  Student,  Troy  High  School,  Troy,  Ml,  June  17  -  September  30,  2013. 

Micro-CT  and  DTS  image  processing  and  analysis. 

16.  Ellen  Yang,  Student,  Troy  High  School,  Troy,  Ml,  July  17  -  September  30,  2013. 

Micro-CT  and  DTS  image  processing  and  analysis. 

17.  Matthew  Varga,  Student,  Wayne  State  University,  BMS  Masters  Program,  June  2013  -  present 
Image  processing  and  analysis. 

18.  Omaima  Bokhari,  Student,  Wayne  State  University,  School  of  Medicine,  Year  4,  August  2013  -  present. 
DTS  image  processing  and  analysis. 

CONCLUSION 

DTS  generated  images  analyzed  by  fractal,  MIL  and  LFD  methods  provided  good  to  excellent  predictions  of  in 
vitro  vertebral  trabecular  bone  stiffness  and  whole  vertebral  strength.  This  technique  could  augment  the  current 
clinical  imaging  modalities  (DEXA  and  CT)  to  increase  the  prediction  accuracy  of  the  fracture  risk  in  patients. 

Future  work  is  planned  to  adapt  the  tomosynthesis  reconstruction  techniques  to  multipurpose  imaging  systems 
so  as  to  achieve  multidirectional  DTS  scan  capabilities  in  vivo  and  extend  the  feasibility  studies  to  clinical 
experiments  where  the  ability  of  DTS  to  predict  fracture  risk  will  be  examined  using  clinical  cohorts  known  to  be 
at  high  risk,  such  as  those  with  prevalent  fractures  or  metabolic  bone  disease.  Comparison  of  DTS  with  other 
clinically  available  modalities  such  as  DEXA  and  CT  are  also  planned. 
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Digital  Tomosynthesis-Derived  Microstructural 
Parameters  Predict  Cancellous  Bone  Stiffness  in 
Human  Vertebrae 

Kim,  Woong1;  Oravec,  Daniel1;  Sander,  Edward2;  Divine,  George 
W1;  Flynn,  Michael  J1;  Yeni,  Yener  N1 


1.  Bone  and  Joint  Center,  Henry  Ford  Hospital,  Detroit,  MI,  United  States. 

2.  Dept  of  Biomedical  Engineering,  The  University  of  Iowa,  Iowa  City,  LA,  United  States. 

Introduction:  Clinical  modalities  currently  available  for  the 
assessment  of  fracture  risk  in  the  spine  largely  rely  on  the  bone 
mineral  density  (BMD)  measurements;  however  their  ability  to 
provide  information  on  the  trabecular  structure  is  limited.  Digital 
tomosynthesis  (DTS)  is  a  relatively  new  imaging  modality  that  has 
many  advantages  over  the  conventional  CT  and  radiography, 
where  DTS  can  generate  multiple  image  slices  with  superior  in¬ 
plane  resolution  than  that  of  CT  but  with  only  1 /6th  of  the  radiation. 
Nekkanty  et  al.  have  recently  reported  that  fractal  dimension  (FD) 
as  measured  from  DTS  images  of  human  vertebrae  is  correlated  to 
stereologic  measures  of  cancellous  microstructure  obtained  from 
p-CT  images  [1].  Also  they  have  shown  that  FD  could  be  used  to 
predict  the  compressive  stress  state  in  the  micro  structure.  The 
objectives  of  this  study  was  to  investigate  if  the  MIL  (mean 
intercept  length)  and  LFD  (line  fraction  deviation  [2])  parameters 
of  DTS  images  obtained  from  human  vertebrae  could  i)  be 
correlated  to  the  p-CT  stereological  parameters  and  ii)  predict 
trabecular  bone  stiffness,  independent  from  bone  volume  fraction. 

Methods:  T6,  T8,  Til,  and  L3  vertebrae  from  5  female  and  5 
male  cadavers  (age  63-90yr)  were  p-CT  scanned  and  reconstructed 
at  45  pm  voxel  size.  In  each  trabecular  bone  region,  mean  (Av) 
and  standard  deviation  (SD)  of  the  bone  volume  fraction  (BV/TV), 
trabecular  thickness  (Tb.Th),  number  (Tb.N),  and  separation 
(Tb.Sp),  and  degree  of  anisotropy  (DA)  were  calculated. 
Cancellous  bone  apparent  modulus  (E)  was  calculated  from  a 
finite  element  (FE)  simulation  of  inferior- superior  compression  [3]. 
The  same  bones  were  scanned  using  DTS  (Shimadzu  Sonialvision 
Safire  II)  and  reconstructed.  In  each  image  set,  the  coronal  central 
slice  was  used  for  MIL  and  LFD  analyses  (Fig  1).  Mean  (Av), 
standard  deviation  (SD),  maximum  (Max)  and  degree  of 
anisotropy  (DA)  were  recorded  (MIL. DA  =  ellipse  fitted  Max/Min 
and  LFD. DA  =  Raw  Max/LFD  orthogonal  to  Max).  Mixed  model 
regressions  of  p-CT  parameters  to  DTS  parameters  were  examined 
(Table  1).  Multi  regression  models  were  constructed  to  examine 
the  relationship  of  MIL  and  LFD  parameters  with  E  with  and 
without  BV/TV.  In  all  regression  models,  the  donor  variable  was 
included  as  a  random  effect  to  account  for  pseudo-replication. 


- MIL  270  - LFD  270 


. Ellipse  fitted  Max  &Mn  . .  Max  &  Qthogonal  LFD 

Figure  1  -  [A]  Coronal  view  of  DTS  image  of  a  L3  vertebra  and 

[B]  region  analyzed  with  [C]  MIL  (mean  intercept  length)  and  [D] 
LFD  (line  fraction  deviation)  methods.  The  overall  alignment  of 
MIL  and  LFD  plots  corresponds  to  that  of  the  trabecular  texture. 
Dotted  red  lines  mark  vectors  of  ellipse  fitted  max  and  min  MIL  in 

[C] ,  and  max  LFD  and  orthogonal  to  the  maximum  LFD  in  [D]. 


Results:  MIL.Max,  MIL.DA,  MIL.SD,  LFD. Max  and  LFD.Av 
from  DTS  were  significantly  related  with  Tb.Th. Av,  BV/TV. SD, 
Tb.Th. SD,  Tb.N.SD  and  DA  from  the  p-CT  stereological 
parameters  (Table  1,  p  <  0.001  to  p  <  0.05).  DTS  parameters  were 
significant  in  multiple  regression  models  of  E  after  accounting  for 
BV/TV  (Fig  2).  Multiple  regression  models  of  E  using  DTS 
parameters  alone  (MIL.DA  +  MIL.Av  +  LFD. Max)  were  as 
predictive  as  BV/TV  alone  (R2adj  =  0.832  vs.  0.819  respectively). 


R 

(p-value) 

MLMax 

ML  DA 

MLSD 

LFD.Max 

LFD.Av 

Tb.Th.Av 

+0.82 

(0.01) 

0.83 

(0.01) 

0.77 

(0.04) 

0.77 

(0.02) 

BV/TV.SD 

-0.86 

(0.02) 

7b.Th.SD 

+0.87 

(0.001) 

0.86 

(0.001) 

0.81 

(0.03) 

7b.N.SD 

-0.86 

(0.001) 

-0.86 

(0.01) 

-0.86 

(0.001) 

DA 

+0.75 

(0.02) 

+0.73 

(0.05) 

Table  1-  R  of  mixed  linear  regression  of  p-CT  (rows)  against  DTS 
parameters  (columns).  Only  those  with  p  <  0.05  are  shown. 
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Figure  2  -  [A]  Actual  vs.  predicted  stiffness  plot  of  the  multi 
regression  model  based  on  both  DTS  and  BV/TV  parameters.  [B, 
C,  D]  The  regression  model’s  leverage  plots  for  the  individual 
effects  (MIL.DA,  LFD.Max  and  BV/TV.Av)  with  the  associated 
p-values.  The  results  show  that  stiffness  prediction  with  the  model 
is  excellent  (R2acij  =  0.947). 

Discussion:  We  demonstrated  that  correlates  of  the  standard  p-CT 
stereological  parameters  can  be  obtained  from  DTS  images  using 
MIL  and  LFD  methods.  Also  we  were  able  to  demonstrate  that 
DTS  images  can  be  used  to  predict  stiffness  of  human  vertebral 
bone.  In  fact,  DTS  derived  model  performed  equally  well  as  the 
BV/TV-only  model  indicating  that  DTS  alone  is  potentially  as 
powerful  a  tool  as  those  measuring  BMD  for  predicting  vertebral 
stiffness.  When  BV/TV.Av  (a  surrogate  of  bone  mass)  and  DTS 
parameters  were  combined,  the  explained  variability  in  E 
increased  from  82%  to  95%,  suggesting  that  DTS  contains 
trabecular  structural  information  that  is  independent  of  BV/TV.Av 
and  therefore  it  could  be  used  as  a  complementary  modality  to 
DEXA  and  tomographic  BMD  measures.  This  in  turn  opens  an 
exciting  possibility  that  DTS  can  be  used,  not  only  as  a  qualitative 
imaging  tool,  but  as  a  quantitative  tool  to  assess  bone  fracture  risk 
while  subjecting  patients  with  a  relatively  low  radiation.  Because 
cancellous  bone  stiffness  is  a  strong  correlate  of  bone  strength,  we 
expect  that  these  methods  will  perform  equally  well  for  predicting 
bone  strength.  However,  it  remains  to  be  determined  to  what 
extent  DTS  derived  parameters  can  predict  whole  bone  strength 
and  fracture  risk  in  vivo. 

Significance:  DTS  generated  images  analyzed  by  MIL  and  LFD 
methods  provided  good  to  excellent  predictions  of  in  vitro 
vertebral  trabecular  bone  stiffness  and  thus  this  technique  could 
augment  the  current  clinical  imaging  modalities  (DXA  and  CT)  to 
increase  the  prediction  accuracy  of  the  fracture  risk  in  patients. 
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References:  [1]  Nekkanty  et  al.,  2011,  57th  ORS,  p.  670.  [2] 
Geraets,  W.  1998,  Bone,  23(4):  p.  383-388.  [3]  Yeni  et  al,  2011, 
Bone,  49:  886-94. 


APPENDIX  2 


2013  Siemens  Competition 

Math  :  Science  :  Technology 

Supplemental  Form 

(Please  save  to  your  hard  drive,  type  in  answers  and  print.  Make  three  copies  and  attach  one  to  each  copy  of  your  Research  Project.) 

Student  Information 


Name  of  Individual  or  Team  Leader: 

Angela  Xiao 

Names  of  Parents  of  Individual  or  Team  Leader: 

Steve  Xiao,  Kathy  Jin 

Name  of  Individual's  or  Team  Leader's  High 

School: 

Troy  High  School 

High  School  City  and  State: 

Troy,  Ml 

Name  of  Team  Member  2: 

Ellen  Yang 

Names  of  Parents  of  Team  Member  2: 

George  Yang,  Mei  Li 

Name  of  Team  Member  3: 

Names  of  Parents  of  Team  Member  3: 

Research  Project  Title: 

Topographical  Measurements  of  Vertebral  Endplate  to  Predict  Bone  Fragility:  Comparison  of  Digital 
Tomosynthesis  and  High  Resolution  Computed  Tomography  with  Microcomputed  Tomography 

Project  Type: 

□  individual  1X1  Two-person  team  1  1  Three-person  team 

Publications 


Has  this  research  project  been  published,  accepted, 
or  submitted  for  publication? 

□  YES  -  published  □  YES  -  accepted  □  YES  -  submitted  1X1  NO 

If  YES,  please  name  the  journal(s) 

Acknowledgements 

You  must  acknowledge  all  individuals  who  assisted  or  advised  you  in  any  way  (i.e.  mentor,  teacher,  family  members,  lab  assistants,  other  students). 
Please  include  the  person's  name,  title,  and  institution  and  describe  their  participation/contribution.  You  are  required  to  complete  the 
Acknowledgements  section  and  must  include  the  name,  title,  institution,  and  role  of  any  individuals  who  assisted  with  the  Research  Project, 
including  family  members.  Receiving  assistance  from  family  members  does  not  affect  the  judging  in  any  way,  but  you  must  acknowledge  their 
support.  Failure  to  identify  others  who  assisted  you  with  any  aspect  of  your  research  may  result  in  disqualification. 


Name,  Title,  Institution,  Participation/Contribution 

SAMPLE  -  Dr.  John  Smith ,  Professor  of  Biology,  Sample  University ,  Assisted  with  the  collection  of  data 


Dr.  Yener  Yeni,  Head  of  Biomechanics  Section,  Henry  Ford  Hospital,  Provided  project  oulines 

Daniel  Oravec,  Research  Engineer,  Henry  Ford  Hospital,  Assisted  with  collection  of  data,  students  were  provided  background 
information  and  extensive  training  on  the  use  of  ImageJ  and  application  of  established  protocols  for  endplate  topography 
measurement.  Provided  outline  for  interpretation  of  statistical  results  and  contributed  toward  discussion. 


Additionally,  please  identify  the  institutions  (i.e.  laboratory,  university,  etc.)  where  you  performed  your  research,  if  not  already  acknowledged 
above. 


Institution 

City 

State 

Henry  Ford  Health  System 

Detroit 

Ml 

May  2013 


Please  expand  to  additional  pages  as  needed. 


Topographical  Measurements  of  Vertebral  Endplate  to  Predict  Bone  Fragility: 
Comparison  of  Digital  Tomosynthesis  and  High  Resolution  Computed 
Tomography  with  Microcomputed  Tomography 

ABSTRACT 

Bone  mineral  density  measurements  help  determine  strength  of  human  vertebral  bodies  and  their 
fracture  potential,  but  BMD  alone  is  not  the  most  accurate  predictor.  Research  indicates  that 
surface  topography  measurements  of  vertebral  endplates  can  increase  the  accuracy  in  assessment 
of  bone  fragility.  Thus,  to  maximize  prediction  accuracy,  clinical  measurements  with  little 
radiation  that  include  superior  and  inferior  endplates  would  be  useful.  Our  objective  was  to 
investigate  the  feasibility  of  measuring  endplate  topography  using  clinically  available  imaging 
methods:  High  resolution  computed  tomography  (HRCT)  and  digital  tomosynthesis  (DTS),  a 
scan  with  a  small  rotation  angle  that  yields  images  of  limited  field  depth,  were  examined  in 
comparison  with  micro  computed  tomography  (pCT),  the  current  gold  standard  for  scanning 
non-living  specimens.  Topographical  endplate  measurements  were  recorded;  five  statistical 
parameters  (average,  standard  deviation,  coefficient  of  variation,  kurtosis,  and  skewness)  were 
calculated  for  each  set  of  scans.  After  statistical  analysis,  a  significant  correlation  was  found  in 
all  parameters  between  DTS  and  pCT,  and  HRCT  and  pCT  for  the  superior  endplate.  This 
demonstrates  that  DTS,  a  scan  with  double  the  resolution  of  CT  scans  in  the  plane  of  endplate 
measurement  and  less  radiation  exposure,  is  a  potentially  viable  approach  to  quantify  endplate 
topography  in  a  clinical  setting. 


EXECUTIVE  SUMMARY 


Osteoporosis  is  a  progressive  disease  where  bone  becomes  more  porous  and  less  dense,  making 
it  weaker  and  more  prone  to  breakage.  It  affects  9  million  Americans  putting  them  all  at  a  risk  of 
severe  bone  fracture.  There  are  two  major  types  of  bone:  cancellous  bone  is  the  inner,  softer 
portion  of  the  bone,  and  cortical  bone  is  the  outer,  harder  layer  of  bone.  Areas  with  a  higher  ratio 
of  cancellous  bone  to  cortical  bone  such  as  the  spine  are  closely  associated  with  increased 
fracture  risk,  because  cancellous  bone  undergoes  turnover  (removal  and  replacement  of  old 
bone)  at  a  faster  rate  than  cortical  bone  does.  Vertebral  endplates  are  composed  of  a  layer  of 
thickened  cancellous  bone,  so  increased  osteoclast  activity  (removal  of  bone  tissue)  can  result  in 
endplate  abnormality  and  dysfunction.  Vertebral  fractures  can  be  detected  based  on  the  presence 
of  endplate  deformities  and  a  lack  of  parallelism  of  endplates.  Thus,  to  increase  accuracy  of 
diagnosis  and  reduce  risk  of  fracture,  scans  should  include  superior  and  inferior  portions  of 
vertebral  endplate.  We  investigated  three  types  of  scans:  microcomputed  tomography  (pCT), 
high  resolution  computed  tomography  (HRCT),  and  digital  tomosynthesis  (DTS);  of  the  three, 
DTS  administers  the  smallest  radiation  dose  (1/5  that  of  HRCT).  pCT  has  the  highest  resolution, 
but  is  not  feasible  for  clinical  use  due  to  high  radiation,  so  it  is  used  only  as  a  standard  of 
comparison  for  the  other  two  scans,  both  of  which  have  been  approved  for  clinical  use.  In  our 
study,  a  set  of  1 17  vertebral  bones  were  scanned  using  each  of  the  three  scanning  mechanisms. 
Endplate  topography  measurements  were  statistically  analyzed  for  each  set  and  five  parameters 
were  calculated  (average,  standard  deviation,  coefficient  of  variation,  kurtosis,  and  skewness).  A 
correlation  was  discovered  between  both  DTS  and  HRCT  in  relation  to  pCT  for  the  superior 
endplate.  Because  of  this,  we  can  infer  that  DTS  is  a  potential  alternative  in  predicting  bone 
fracture,  because  of  its  accuracy  in  endplate  topography  measurement  and  its  low  radiation  dose. 


Topographical  Measurements  of  Vertebral  Endplate  to  Predict  Bone 
Fragility:  Comparison  of  Digital  Tomosynthesis  and  High  Resolution 

Computed  Tomography  with  Microcomputed  Tomography 

I.  INTRODUCTION 

Osteoporosis  is  a  bone  disease  that  progresses  silently,  often  with  few  symptoms;  as  time  goes 
on,  the  risk  of  fracture  significantly  increases.  Identifying  and  treating  patients  at  risk  of  fracture 
can  substantially  reduce  the  long  term  burden  of  osteoporosis,  reducing  the  risk  of  first  fracture 
from  8%  to  2%  can  reduce  the  5-year  fracture  incidence  from  approximately  34%  to  10% 
(Lindsay  et  al,  2005).  Bone  mineral  density  (BMD)  is  commonly  used  in  assessment  of  bone 
quality  and  fracture  risk  (Cranney  et  al.  2007),  but  taken  alone,  BMD  may  not  be  the  best  single 
determinant  of  whole  bone  mechanical  properties  (Cheng  et  al.,  1997;  Singer  et  al.,  1995). 

Recent  studies  have  shown  that  topographical  features  of  human  vertebral  body  endplates  are 
associated  with  vertebral  bone  mechanical  properties  and  may  be  important  determinants  of 
bone  fragility  (Eswaran  et  al.,  2007;  Langrana  et  al.,  2006;  Nekkanty  et  al.,  2010).  Nekkanty  et 
al.  found  that  kurtosis  of  the  endplate  topography  statistical  distribution,  which  captures 
topographical  features  such  as  peaks  or  deep  valleys,  are  significantly  correlated  at  the  inferior 
endplate  with  finite  element-derived  mechanical  parameters  such  as  work  to  fracture, 
displacement  at  maximum  load,  and  strain  to  maximum  load.  In  such  studies  of  bone 
microstructure,  microcomputed  tomography  (pCT)  is  typically  used  as  a  gold  standard  as  it 
provides  a  high  resolution  3 -dimensional  depiction  of  bone  microarchitecture.  However,  pCT  is 
only  suitable  in  a  laboratory  setting  with  in  vitro  specimens  due  to  size  limitations  and  high 
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radiation  dose.  To  date,  the  use  of  clinical  imaging  modalities  in  assessment  of  endplate 
topography  has  not  been  studied. 

Computed  tomography  (CT)  is  used  clinically  in  assessment  of  bone  structure,  but  in-plane 
reconstructed  pixel  size  is  typically  on  the  order  of  0.3-1  mm  with  a  resolution  of  approximately 
1  mm  (Thornton  et  al,  2006).  Recent  advances  in  flat  panel  digital  detectors  have  realized  a 
technology  called  digital  tomosynthesis  (DTS),  a  scan  that  administers  a  significantly  lower 
dose  of  radiation  while  maintaining  high  in  plane  resolution.  Digital  tomosynthesis  is  a 
tomographic  imaging  modality  in  which  a  series  of  projection  images  are  acquired  over  a 
limited  arc,  with  the  x-ray  source  pivoting  and  translating  opposite  the  direction  of  a  flat  panel 
detector  encased  in  the  scanning  bed.  DTS  delivers  approximately  1/5  the  dose  of  CT  and 
provides  good  resolution  of  cancellous  texture  and  cortical  bone  geometry  (Kim  et  al.,  2013, 
Nekkanty  et  al.,  2011).  Tomosynthesis  reconstructions  are  formed  in  the  lateral-medial  or 
anterior-posterior  direction,  in  contrast  with  CT  which  produces  an  axial  image  with  high  slice 
thickness  along  the  superior- inferior  direction.  Topographic  features  of  the  endplate  surface  are 
thus  captured  in-plane  with  a  pixel  size  of  approximately  0.3  mm  in  DTS,  rather  than  out-of- 
plane  with  0.75  mm  slice  thickness  in  an  axial  HRCT  image. 

In  this  study,  endplate  topography  measurements  of  digital  tomosynthesis  (DTS)  and  high 
resolution  computed  tomography  (HRCT)  were  compared  to  similar  measurements  from  pCT  in 
order  to  assess  the  viability  of  endplate  topography  measurement  in  a  clinical  imaging  modality. 
Due  to  better  in-plane  resolution  of  the  vertebral  endplate,  we  hypothesized  that  endplate 
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characteristics  measured  from  tomosynthesis  images  would  be  better  correlated  with  pCT  than 
measurements  from  HRCT  (Figure  1). 


Figure  1:  Comparison  of  images  taken  from  sagittal  plane  in  similar  regions  from  (a)  DTS,  (b) 
HRCT,  and  (c)  pCT.  Images  are  resized  to  show  detail  at  comparable  scale. 

II.  METHODS 
2.1  Subjects 

Human  cadaveric  thoraco-lumbar  spines  were  acquired  under  local  IRB  approval  from  tissue 
banks,  and  four  vertebral  bodies  (T10,  T1 1,  T12,  LI)  were  harvested  from  42  donors  of  varying 
gender  and  age.  Donors  with  HIV,  hepatitis  or  infectious  diseases,  cause  of  death  involving 
trauma  or  a  history  of  spine  surgery,  glucocorticoid  use,  or  metabolic  diseases  known  to  affect 
bone  were  not  included.  Vertebral  bodies  were  dissected,  soft  tissue  and  posterior  elements  were 
removed,  and  specimens  were  stored  wrapped  in  saline-soaked  gauze  at  -20°  C  until  scanning 
was  performed.  The  donor  set  consisted  of  23  men  and  19  women  all  between  the  ages  of  36  and 
100.  Collectively,  the  vertebrae  of  these  donors  formed  a  set  of  1 17  bones. 
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2.2  HRCT,  DTS,  and /uCT  Scanning  Procedures 


Specimens  were  mounted  and  consistently  aligned  in  a  custom  radiolucent  scanning  tank  filled 
with  0.9%  saline  and  scanned  using  DTS  (Shimadzu  Sonialvision  Safire  II)  and  high  resolution 
CT  (Siemens  Sensation  64)  (Figure  2).  Alignment  was  ensured  using  a  radiolucent  clamping 
system  such  that  the  anterior-posterior  and  lateral-medial  anatomical  directions  corresponded  to 
the  reconstruction  x  and  y  axes,  respectively.  The  same  specimens  were  scanned  using  a  custom- 
built  pCT  system  and  reconstructed  at  an  isotropic  voxel  size  of  40  micrometers.  The  pCT 
system  used  in  the  study  was  home  built  based  on  the  hardware,  data  acquisition,  and 
reconstruction  methods  of  the  pCT  system  have  been  previously  described  (Reimann  et  al., 
1997).  The  presently  operating  system  uses  a  Kevex  16-watt  x-ray  source  with  a  9-micron  focal 
spot,  a  1888xl408-pixel  Varian  PaxScan  2520  flat  panel  digital  x-ray  detector  with  127-micron 
pitch,  a  Newport  precision  rotational  stage,  and  control  software  running  under  Windows  XP. 


(a)  (b)  (c) 


Figure  2:  Vertebral  bodies  were  mounted  to  trays  via  a  precise  clamping  system  contained 
within  radiolucent  boxes  fdled  with  saline.  Boxes  were  scanned  in  digital  tomosythesis  (a)  and 
computed  tomography  systems  (b).  The  same  vertebral  bodies  were  positioned  in  a  radiolucent 
cylinder  and  scanned  using  microcomputed  tomography  (c). 
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2.3  HRCT,  DTS,  and /uCT Data  Collection 


For  HRCT  images,  a  single  threshold  value  in  Hounsfield  Units  was  manually  determined  as  the 
minimum  value  that  delineates  bone  from  soft  tissue.  The  threshold  value  was  used  in  a  custom 
segmentation  algorithm  within  an  in-house  image  processing  suite  which  produced  a  binarized 
volume  (Zauel  et  al.,  2006).  After  thresholding,  the  disconnected  voxels  within  the  volume  were 
removed  via  a  continuity  check  algorithm.  Volume  masks  separating  cortical  and  cancellous 
bone  were  segmented  using  a  previously-described  semi-automatic  method  (Oravec  et  al.,  2012, 
Buie  et  al.  2007).  The  segmentation  algorithm  consists  of  dilating  the  binarized  vertebral  image 
twice  (closing  porosity  within  and  on  the  surface  of  the  vertebra),  applying  a  median  filter 
(connecting  surfaces  and  smoothing  processing  artifacts),  and  eroding  back  twice.  The  resulting 
volume,  a  solid  mask  image  of  the  whole  vertebral  body,  was  further  cropped  into  separate 
volumes  representing  superior  and  inferior  endplates.  TopoJ  (Hovis,  2013),  a  plugin  written  for 
ImageJ  (Rasband,  1997-2013),  was  used  to  analyze  the  topography  of  superior  and  inferior 
endplates  separately.  TopoJ  calculates  a  2D  height  map,  with  each  pixel  in  the  map  representing 
the  depth  from  a  fixed  plane  (the  first  slice)  of  the  first  encounter  of  bone  along  the  superior- 
inferior  axis.  The  depth  distribution  (background  and  holes  were  eliminated  from  the  analysis) 
was  recorded  into  a  text  file  in  which  each  row  represented  a  single  pixel  in  2D  height  map. 
Distribution  statistics  such  as  average  (AV),  standard  deviation  (SD),  coefficient  of  variation 
(CV),  skewness  (Skew),  and  kurtosis  (Kurt)  were  calculated  from  HRCT  depth  measurements. 

A  similar  process  was  performed  for  DTS  images.  A  global  threshold  value  was  manually 
determined  to  delineate  bone  from  soft  tissue  and  air.  Due  to  a  blurring  effect  in  the  highest  and 
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lowest  slices  of  the  DTS  reconstructions,  a  central  substack  of  25  slices  (25  mm)  was  created  and 
the  image  was  binarized  in  ImageJ  using  the  recorded  threshold  value.  Binarized  images  were 
cropped  into  separate  volumes  representing  superior  and  inferior  endplates.  TopoJ  was  again 
used  to  calculate  depth  distributions  for  DTS  endplate  images.  The  same  distribution  statistics 
were  calculated  from  DTS  depth  measurements. 

Threshold  values  for  pCT  reconstructions  were  calculated  from  a  300x300x300  voxel  central 
cube  of  cancellous  bone  using  the  Otsu  method  in  ImageJ  and  the  image  was  binarized  using  this 
global  threshold  value.  A  custom  algorithm  was  used  to  produce  a  depth  frequency  distribution 
representing  distances  from  the  superior-most  and  inferior-most  reference  planes  to  the  inner 
surfaces  of  the  vertebral  endplate.  The  same  distribution  statistics  were  calculated  using  the 
depth  frequency  distribution  produced  by  this  method  for  pCT. 

2.3  Statistical  analysis 

The  main  interest  was  in  the  relationship  between  similar  measurements  from  the  DCT  and  DTS 
images  and  between  those  from  the  pCT  and  high  resolution  CT  (HRCT)  images  as  determined 
using  regression  models.  The  analysis  started  with  examination  of  interactions  between  vertebral 
level  and  the  predictor  of  interest.  Mixed  models  were  constructed  using  one  of  pCT  variables  as 
the  outcome,  the  corresponding  DTS  or  HRCT  variable,  vertebral  level  and  the  interaction  of  the 
latter  two  as  the  effect  variables.  The  model  also  included  a  random  subject  variable  to  account 
for  pseudoreplication  due  to  using  multiple  vertebral  levels  from  some  subjects.  If  a  significant 
interaction  was  found,  the  data  were  split  by  vertebral  level  and  a  separate  simple  regression 
analysis  was  performed  between  the  pCT  and  the  DTS  (or  HRCT)  variable.  Nonsignificant 


6 


interaction  terms  were  removed  and  the  models  were  rerun  with  the  main  effects  until  no 


nonsignificant  effect  was  present.  If  the  data  could  be  pooled  through  this  elimination  procedure 
and  a  significant  positive  relationship  remained  between  the  endplate  variables  from  the  two 
methods,  such  measurements  were  deemed  most  successful. 

Paired  t-tests  were  also  run  to  assess  differences  in  mean  of  statistical  distribution  parameters 
between  superior  and  inferior  endplates.  All  analyses  were  performed  in  JMP  (SAS  Institute  Inc., 
Cary,  NC). 

III.  RESULTS 

Statistical  Results  are  summarized  in  Table  1.  Mixed  models  indicated  significant  interactions 
between  pCT  and  DTS  for  superior  SD  (p<0.03)  and  skewness  (p<0.02),  as  well  as  between  pCT 
and  HRCT  for  inferior  AV  (p<0.05).  When  vertebrae  were  separated  by  level,  significant 
positive  correlations  were  found  for  DTS  superior  SD  and  skewness  at  all  four  levels  (T10,  Til, 
T12,  and  LI),  so  the  levels  were  pooled  in  a  unified  mixed  model  in  which  both  SD  and 
skewness  were  found  to  be  significant  (p<0.0001,  R2Adj=0.21;  p<0.0001,  R2Adj=0.53).  When 
separated  by  level,  HRCT  average  for  inferior  endplate  was  significant  only  for  T10-T12. 

At  the  inferior  endplate,  only  significant  level  effects  were  seen  in  DTS  skewness  (p<0.04)  and 
HRCT  kurtosis  (p<0.03)  and  skewness  (p<0.04).  Neither  level  nor  variable  were  significant  for 
models  of  superior  DTS  CV  (p>0.3)  and  inferior  CT  SD  (p>0.3)  and  CV  (p>0.08). 
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In  the  case  that  both  variable  and  level  were  found  to  be  significant,  mixed  models  were  rerun 
removing  level  as  an  effect.  Superior  DTS  CV  was  significant  without  level  (p<0.0001,  R  Adj 
=0.17),  as  was  HRCT  SD  (p<0.0001,  R2Adj  =0.19)  and  inferior  DTS  kurtosis  (p<0.02,  R2Adj 
=0.03).  Variables  which  remained  significant  following  the  elimination  procedure  (no  significant 
interaction  or  level  effects)  were  superior  DTS  AV  (pO.OOOl,  R  Adj=0.21)  and  kurtosis 
(pO.OOOl,  R2Adj  =0.25),  superior  HRCT  AV  (pO.OOOl,  R2Adj  =0.37),  CV  (p<0.0001,  R2Adj 
=0.41),  skewness  (p<0.002,  R2Adj  =0.58)  and  kurtosis  (p<0.0007,  R2Adj  =0.41),  and  inferior  DTS 
AV  (pO.OOOl,  R2Adj  =0.46)  and  SD  (p<0.02,  R2Adj  =0.29).  There  were  no  significant  variables  at 
the  inferior  endplate  for  CT  without  separation  by  level. 

Paired  t-test  results  indicated  that  the  mean  values  of  average,  standard  deviation,  coefficient  of 
variation  and  kurtosis  of  endplate  topography  distributions  measured  from  pCT  images  are 
significantly  different  between  superior  and  inferior  endplates.  Average  and  Kurtosis  were 
significantly  higher  at  the  superior  endplate,  while  standard  deviation  and  coefficient  of  variation 
were  lower. 
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Table  1:  Statistical  results  comparing  pCT  to  DTS  and  HRCT  topography  measurements  at 
superior  and  inferior  endplates.  All  significant  parameters  had  positive  correlations.  If  a  model 
was  rerun  with  a  non-significant  interaction  term  removed  and  both  level  and  variable  were 
significant,  the  model  was  rerun  without  level  (signified  as  *  Variable).  In  the  case  of  a 
significant  interaction,  if  all  4  levels  were  separated  by  level  and  found  to  be  significant  and 
positive,  levels  were  pooled  in  a  unified  mixed  model  (signified  as  **Variable). 


Comparison 

Result 

Y  vs 

X 

mCTAV 

DTSAV 

p<0.0001,R2Adj=0.21 

mCTSD 

DTSSD 

**p<0.0001,  R2Adj  =0.27 

mCTCV 

DTSCV 

*  p<0.0001,  R2Adj  =0.17 

mCTSkew 

DTS_Skew 

**  p<0.0001,  R2Adj  =0.53 

_o 

<D 

mCTKurt 

DTSKurt 

pO.OOOl,  R2Adj  =0.25 

GO 

mCTAV 

HRCTAV 

p<0.0001,  R2Adj  =0.37 

mCTSD 

HRCT  SD 

*  p<0.0001,  R2Adj  =0.19 

mCTCV 

HRCTCV 

p<0.0001,  R2Adj  =0.41 

mCT_Skew 

HRCT  Skew 

p=0.0014,  R2Adj  =0.58 

mCTKurt 

HRCTKurt 

p=0.0006,  R2Adj  =0.41 

mCTAV 

DTSAV 

pO.OOOl,  R2Adj  =0.46) 

mCTSD 

DTSSD 

p=0.0136,  R2Adj  =0.29) 

mCTCV 

DTSCV 

N.S. 

mCTSkew 

DTSSkew 

N.S. 

_o 

mCTKurt 

DTSKurt 

*p=0.0174,  R2Adj  =0.03 

?— 1 

G-h 

a 

mCTAV 

HRCTAV 

T10  (p=0.0034,  R2Adj  =0.23),  Til  (pO.OOOl, 

i— i 

R2Adj  =0.23),  T12  (p=0.0297,  R2Adj  =0.15) 

mCTSD 

HRCTSD 

N.S. 

mCTCV 

HRCTCV 

N.S. 

mCTSkew 

HRCTSkew 

N.S. 

mCTKurt 

HRCTKurt 

N.S. 
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Table  2:  Paired  t-test  results  comparing  statistical  distribution  parameters  between  superior  and 
inferior  endplates  in  pCT  images.  Average  parameters  for  superior  and  inferior  endplate  are 
presented  as  ±  standard  deviation. 


Parameter 

p-value 

Superior  Endplate 

Inferior  Endplate 

AV  (mm) 

<0.0001 

3.832  ±  1.367 

2.956  ±  1.148 

SD  (mm) 

<0.0001 

1.205  ±0.523 

1.891  ±0.908 

CV 

<0.0001 

0.323  ±  0.086 

0.734  ±0.535 

Skew 

0.4358 

-0.598  ±0.589 

-0.664  ±  1.022 

Kurt 

<0.0001 

3.428  ±  1.288 

2.226  ±  0.898 

IV.  DISCUSSION 

To  our  knowledge,  this  is  the  first  study  to  demonstrate  the  ability  to  measure  vertebral  endplate 
topography  using  clinical  imaging  modalities.  The  authors  successfully  prepared  DTS  and  HRCT 
images  for  analysis,  performed  voxel-based  endplate  topography  measurements,  and  analyzed 
the  topographical  distribution  for  comparison  with  similar  measurements  taken  from  pCT. 

The  level  of  success  was  promising,  in  that  both  tested  modalities  showed  statistically  significant 
correlations  with  pCT  approaching  a  moderate  level  of  agreement  at  the  superior  endplate  for  all 
measured  parameters  (R2Adj=0. 17-0.58).  However,  at  the  inferior  endplate,  average  and  standard 
deviation  were  moderately  correlated  while  higher  order  statistical  moments  such  as  skewness 
and  kurtosis  were  not  significantly  correlated  with  pCT. 
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The  lack  of  ability  to  measure  higher  order  statistical  distribution  parameters  at  the  inferior 
endplate  may  be  explained  as  a  function  of  measurement  resolution.  Kurtosis  of  a  topography 
distribution  gives  a  measure  of  its  sharpness  and  is  defined  by  the  presence  of  outliers.  For 
example,  analysis  of  a  relatively  flat  endplate  with  a  small  number  of  deep  valleys  would  yield  a 
low  kurtosis,  while  an  endplate  with  a  high  sloped  rim  would  yield  a  higher  kurtosis.  Likewise, 
characteristic  features  such  as  a  high  remnant  of  a  posterior  element  or  a  bony  spur  are  captured 
by  higher  order  statistical  distribution  parameters.  If  measurement  resolution  is  poor  and  unable 
to  detect  such  outliers,  measurement  of  values  such  as  kurtosis  and  skewness  would  be  grossly 
affected. 

Paired  t-tests  indicated  that  average  depth  from  the  reference  plane  was  significantly  greater  at 
the  superior  endplate.  Provided  that  DTS  has  more  than  seven  times  lower  measurement 
resolution  than  pCT,  ability  to  measure  sharpness  of  the  topography  distribution  will  be 
potentially  largely  affected  by  average  depth.  With  a  greater  average  depth,  there  are  more  steps 
available  at  a  given  voxel  size  to  define  the  characteristic  features  contributing  to  changes  in 
higher  order  statistical  parameters.  Interestingly,  while  all  measured  parameters  in  both 
modalities  were  significantly  correlated  with  pCT  at  the  superior  endplate,  only  DTS  parameters 
were  significantly  correlated  with  pCT  at  the  inferior  endplate.  This  is  consistent  with  the  fact 
that  CT  has  considerably  poorer  resolution  in  the  measurement  direction,  approximately  half  that 
of  DTS  (superior-inferior  pixel  size  of  0.75  mm  vs  0.3  mm).  The  HRCT  images  used  in  this 
study  were  reconstructed  using  a  high  resolution  protocol  and  typical  musculoskeletal  CT  scans 
(slice  thickness  1-3  mm)  would  yield  a  greater  discrepancy  in  results. 


11 


The  use  of  clinical  imaging  modalities,  specifically  DTS,  is  encouraging  and  should  be  explored 
further.  We  intend  to  look  for  options  to  enhance  precision  of  endplate  topography  measurement 
in  tomosynthesis.  As  mentioned  previously,  one  unavoidable  artifact  in  DTS  imaging  is  the 
blurring  of  the  lowest  and  highest  slices  in  the  acquisition,  which  limits  analysis  to  a  central  slab 
of  well-resolved  slices.  In  this  sense,  the  full  vertebral  endplate  was  not  analyzed  in  the  DTS 
image.  Only  lateral-medial  acquisitions  were  assessed  in  this  study.  The  extent  to  which  anterior- 
posterior  acquisitions  can  improve  the  correlations  to  pCT  results  remains  to  be  determined. 

V.  CONCLUSION 

In  this  study,  moderate  correlations  were  found  between  pCT  and  the  two  tested  clinical 
imaging  modalities  (DTS  and  HRCT)  at  the  superior  endplate  for  five  parameters  (average, 
standard  deviation,  coefficient  of  variation,  kurtosis,  and  skewness)  derived  from  endplate 
topography  measurements  .  In  the  statistical  model,  only  DTS  scans  produced  significant  but 
poor  to  moderate  correlations  with  pCT  at  the  inferior  endplate  for  several  parameters.  Lack  of 
correlation  between  pCT  and  HRCT  at  the  inferior  endplate  is  likely  due  to  a  combination  of 
poor  resolution  and  a  lack  of  characteristic  deep  features  which  higher  order  statistics  would 
identify. 

One  of  the  goals  of  this  study  was  to  seek  out  a  scanning  modality  safe  for  clinical  use  (limited 
radiation  dosage),  while  maintaining  high  resolution  images  from  which  endplate  surface 
topography  measurements  can  be  made.  It  was  hypothesized  that  DTS,  with  approximately  1/5 
the  radiation  dose  and  greater  in-plane  resolution  than  CT,  would  be  an  ideal  clinical  imaging 
modality  for  assessing  endplate  topography,  which  has  been  shown  in  pCT  studies  to  correlate 
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with  mechanical  parameters  related  to  bone  fragility.  As  measurements  taken  from  DTS 
produced  poor  correlation  with  comparable  measurements  from  pCT  (considered  the  gold 
standard  for  scanning),  we  conclude  that  further  studies  are  warranted  before  employing  DTS 
scans  in  a  clinical  setting  for  the  purpose  of  endplate  analysis. 

Several  limitations  of  the  present  study  have  been  presented.  Manual  selection  of  threshold 
values  produced  visually  acceptable  segmentation,  however  methods  may  be  refined  to  use 
adaptive  thresholding  techniques  in  order  to  optimize  characterization  of  the  endplate  surface.  In 
addition,  including  anterior-posterior  scans  in  statistical  models  could  provide  missing 
information  due  to  scanning  limitations  in  a  single  direction. 

This  study  demonstrated  the  ability  to  process  DTS  and  HRCT  images  for  the  purpose  of 
performing  endplate  topography  measurements.  Results  are  encouraging  and  point  to  the 
importance  of  continued  investigation. 
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APPENDIX  3 


Digital  Tomosynthesis  of  Human  Vertebral  Bone:  The  Effect  of  Positioning  and  Scan  Orientation  on  Prediction  of  Cancellous  Bone  Stiffness. 

Kim,  Woong;  Oravec,  Daniel;  Divine,  George  W;  Flynn,  Michael  J;  Yeni,  YenerN 
Bone  and  Joint  Center,  Henry  Ford  Hospital,  Detroit,  MI,  United  States. 

Introduction:  Clinical  modalities  currently  available  for  the  assessment  of  fracture  risk  in  the  spine  largely  rely  on  bone  mineral  density  (BMD) 
measurements;  however,  their  ability  to  provide  information  on  the  trabecular  structure  is  limited.  Digital  tomosynthesis  (DTS)  is  a  relatively  new  imaging 
modality  that  has  several  advantages  over  the  conventional  CT  and  radiography,  where  DTS  can  generate  stack  images  of  a  vertebra  with  an  in-plane 
resolution  superior  to  that  of  CT  but  with  only  l/5th  of  the  radiation.  It  has  been  recently  reported  that  fractal  dimension  (FD),  mean  intercept  length  (MIL  ) 
and  line  fraction  deviation  (LFD  [1])  parameters  obtained  from  DTS  images  of  human  vertebrae  are  correlated  to  3D  microstructure  and  stiffness  of 
trabecular  bone  as  determined  from  microcomputed  tomography  (pCT)  [2,  3].  These  correlations  were  achieved  by  using  a  single  image  slice  selected  at  the 
center  of  a  DTS  stack  (20  to  50  slices).  The  scan  geometry  was  such  that  the  movement  of  the  x-ray  source  and  the  detector  were  transverse  to  the  axis  of  the 
vertebra  and  DTS  images  were  acquired  in  an  anteroposterior  (AP)  view  (resulting  in  a  stack  of  coronal  images).  The  objective  of  the  current  study  was  to 
investigate  the  extent  to  which  scan  orientation  (transverse,  parallel  or  oblique  to  the  vertebral  axis)  and  view  (AP  or  lateral)  affect  the  prediction  of 
trabecular  bone  stiffness,  independent  from  bone  volume  fraction.  The  results  are  expected  to  enable  informed  decisions  about  patient  positioning  during 
translation  of  DTS  to  in  vivo  applications. 

Methods:  Thoracic  6,  T8,  Til,  and  L3  vertebrae  from  5  female  and  5  male  cadavers  (age  63-90yr)  were  pCT  scanned  and  reconstructed  at  45  pm  voxel  size. 
In  each  trabecular  bone  region,  mean  bone  volume  fraction  (BV/TV)  was  obtained.  Cancellous  bone  apparent  modulus  (EFe)  was  calculated  from  a  finite 
element  (FE)  simulation  of  inferior-superior  compression  using  heterogeneous  material  moduli  scaled  with  pCT  gray  levels  [2].  The  same  bones  were 
scanned  using  DTS  (Shimadzu  Sonialvision  Safire  II)  in  the  AP  (producing  a  stack  of  coronal  plane  images)  and  lateral  (LM,  producing  a  stack  of  sagittal 
plane  images)  views  while  aligned  axially  (0°),  transversely  (90°)  or  obliquely  (23°)  to  the  superior-inferior  axis  of  the  vertebrae  (a  total  of  6  scans  per  spine). 
The  vertebrae  were  immersed  in  water  during  scans  to  simulate  soft  tissue. 

A  cuboidal  volume  of  interest  was  cropped  out  of  each  reconstructed  DTS  image-stack  to  include  as  much  cancellous  bone  as  possible,  similar  to  the 
volumes  of  interests  (VOIs)  from  the  pCT  images  used  for  the  FE  analysis  (Fig.  1).  Each  cropped  region  was  analyzed  using  the  fractal,  MIL  and  LFD 
methods.  In  fractal  analysis,  fractal  dimension  (FD),  mean  lacunarity  (A)  and  the  slope  of  lacunarity  vs.  box  size  relationship  were  calculated  [4,  5].  In 
MIL  and  LFD  analyses,  the  mean  (Av),  standard  deviation  (SD),  maximum  (Max)  and  degree  of  anisotropy  (DA,  ratio  of  the  principal  measurements  of  MIL 
or  LFD)  were  recorded  for  each  slice  (Fig  2).  The  central  slice  from  each  stack  was  also  analyzed  separately  for  a  comparison  between  the  2D  and  3D 
approaches. 

Stepwise  forward  regression  method  was  used  to  construct  mixed  multiple  regression  models  that  examined  the  relationship  of  fractal,  MIL  and  LFD 
parameters  with  Efe.  In  all  regression  models,  a  subject  variable  was  included  as  a  random  effect  to  account  for  pseudo -replication.  BV/TV  was  introduced 
first  and  forced  to  stay  in  the  model  in  order  to  examine  the  effect  of  DTS  independently  from  bone  mass.  A  separate  analysis  was  performed  for  each  scan 
configuration.  Models  with  high  multicolinearity  (Variance  Inflation  Factor  >  5)  were  rejected. 

Results:  The  models  constructed  using  DTS  parameters  only  (without  BV/TV)  ranged  from  nonsignificant  to  an  explanatory  capability  (R2adj)  of  0.89,  with 
only  the  LM-oblique  (R2adj  =0.85)  and  the  LM-  transverse  (R2adj  =0.89)  cases  reaching  a  value  above  that  is  explained  by  BV/TV  alone  (R2adj  =0.82).  When 
BV/TV  was  present,  all  scan  configurations  provided  parameters  that  increased  the  explained  variability  (R2adj)  in  cancellous  bone  stiffness  over  that  is 
explained  by  BV/TV  (Table  1).  An  exception  is  the  oblique  scans  in  the  lateral  view  in  which  case  the  DTS  parameters  were  either  nonsignificant  (central 
slice)  or  did  not  improve  R2adj  compared  to  BV/TV  alone.  Models  from  transverse  scans  and  the  LM  view  appeared  to  perform  better  than  those  from  the 
other  configurations.  Parameters  related  to  cancellous  bone  anisotropy  (Av(MIL.DA)  or,  similarly,  Av(MIL.SD))  were  persistently  present  in  the  models 
from  axial  and  transverse  scans.  When  improvement  in  R2adj  was  observed  compared  to  the  single-slice  case,  the  model  included  a  term  representing  the 
slice-to-slice  heterogeneity  of  the  microstructure. 

Discussion:  By  using  the  entire  image  stack,  we  were  able  to  calculate  parameters  that  represent  the  3D  heterogeneity  of  the  micro  structure.  Most 
explanatory  models  of  stiffness  had  one  or  more  of  these  heterogeneity  parameters,  demonstrating  a  major  advantage  of  DTS  over  2D  projection  methods.  A 
predictor  set  that  contains  BV/TV,  a  measure  of  anisotropy  and  a  measure  of  heterogeneity  is  consistent  with  results  from  laboratory-scale  high  resolution 
imaging  studies  (pCT)  [6].  Overall,  the  results  are  encouraging  for  further  development  of  bone  assessment  tools  using  DTS. 

The  most  explanatory  model  from  each  scan  configuration  included  different  predictor  parameters.  Different  parts  of  the  cancellous  bone  in  the  analysis 
volume  and  the  dependence  of  DTS  resolution  on  the  orientation  of  structural  features  relative  to  the  scan  orientation  [7]  may  explain  these  differences.  The 
findings  that  DTS  parameters  derived  from  the  LM  views  could  outperform  BV/TV  in  predicting  bone  stiffness  and  that  the  transverse  scans  from  the  LM 
view  resulted  in  more  explanatory  models  in  the  presence  of  BV/TV  suggest  that  the  transverse  orientation  in  the  LM  view  is  the  preferred  configuration  for 
the  best  assessment  of  cancellous  bone  if  a  single  DTS  scan  is  to  be  considered.  However,  it  is  likely  that  scans  in  additional  orientations  or  views  will 
provide  information  complementary  to  that  from  the  LM -transverse  scan  and  increase  prediction  accuracy  for  the  fracture  risk  of  a  whole  vertebra  under  a 
more  complex  loading  than  uniaxial  compression. 

It  is  not  feasible  in  the  current  Safire  system  to  position  a  human  subject  with  a  tilt  more  than  23  with  respect  to  the  scan  table.  However,  a  multi-orientation 
DTS  scan  with  high  resolution  can  be  obtained  (without  sacrificing  the  low -radiation  advantage)  using  a  cone-beam  system  in  a  limited  angle  tomography 
with  appropriate  modifications  to  acquisition  and  reconstruction  protocols.  Current  efforts  are  devoted  to  realization  of  such  modifications. 

Significance:  DTS  generated  image  stacks  analyzed  by  fractal,  MIL  and  LFD  methods  provided  good  to  excellent  predictions  of  in  vitro  vertebral  trabecular 
bone  stiffness,  independent  from  bone  mass.  An  orientation  perpendicular  to  the  axis  of  the  spine  in  an  LM  view  is  the  preferred  configuration  if  a  single 
scan  is  to  be  considered.  Complementary  microstructural  information  can  be  obtained  by  scanning  in  multiple  orientations.  Further  development  for  a 
multidirectional  DTS  is  feasible  and  may  be  necessary. 
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Table  1  -  Summary  of  the  most  explanatory  mixed  multiple  regression  models  (E  =BV/TV  +  DTS)  for  each  scan  and  analysis  combination.  Each  cell  shows 
R2adj  followed  by  the  predictor(s)  in  the  model  (p-value).  Only  those  with  p  <  0.05  are  shown.  BV/TV  is  present  in  all  models  (p<0.0001)  and  therefore  not 
shown  on  the  predictor  list  for  clarity. 
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Figure  1  -  Illustration  of  AP  and  LM  VOIs  in  a  vertebra.  The  red  and  the  blue  lines  represent  the  image  stacks  for  the  LM  and  the  AP  view,  respectively. 
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Figure  2  -  [A]  A  square  region  cropped  out  of  the  DTS  image  is  analyzed  using  the  fractal,  MIL  and  LFD  analysis  methods.  MIL  and  LFD  are  calculated  in 
all  directions  in  the  image.  [B]  A  typical  polar  plot  representation  of  the  LFD  data,  from  which  an  average,  standard  deviation,  max  and  anisotropy  (LFD  at 
90°  /  LFD  at  0°)  can  be  calculated  for  each  slice.  [C]  Each  slice  in  the  stack  is  analyzed.  Stack  average  and  standard  deviation  (e.g.  SD  (LFD.SD))  are 
calculated  as  measures  of  3D  microstructure  within  the  bone. 
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Digital  Tomosynthesis  for  Prediction  of  Human  Whole  Vertebral  Stiffness 

Kim,  Woong;  Oravec,  Daniel;  Maatman,  Thomas;  Divine,  George  W;  Flynn,  Michael  J;  Yeni,  YenerN 
Bone  and  Joint  Center,  Henry  Ford  Hospital,  Detroit,  MI,  United  States. 

Introduction:  Digital  tomosynthesis  (DTS)  has  been  proposed  as  an  imaging  approach  capable  of  providing  information  on  cancellous  bone  micro  structure 
and  thus  as  a  potential  means  for  an  accurate  assessment  of  vertebral  bone  strength  and  fracture  risk  [1,2].  Although  associations  of  DTS  derived  parameters 
with  cancellous  bone  stiffness  have  been  demonstrated,  the  extent  to  which  DTS  parameters  are  associated  with  whole  vertebral  mechanical  properties  is 
unknown.  The  objective  of  the  current  study  was  to  determine  relationships  of  DTS  derived  cancellous  bone  and  vertebral  geometry  properties  with  vertebral 
stiffness.  Large  scale  finite  element  (FE)  estimations  of  vertebral  stiffness  based  on  high  resolution  microcomputed  tomography  (pCT)  images  are  used  as  a 
nondestructive^  determined  surrogate  for  vertebral  strength  due  to  strong  correlations  between  vertebral  stiffness  and  strength. 

Methods:  T6,  T8,  Til,  and  L3  vertebrae  from  5  female  and  5  male  cadavers  (age  63-90yr)  were  p-CT  scanned  and  reconstructed  at  45  pm  voxel  size.  In 
each  trabecular  bone  region,  bone  volume  fraction  (B  V/TV)  was  obtained  and  stiffness  of  whole  bone  (S)  was  calculated  from  a  finite  element  (FE) 
simulation  of  inferior-superior  compression  using  procedures  described  previously  [3,4].  The  same  bones  were  scanned  using  DTS  (Shimadzu  Sonial vision 
Satire  II)  in  anterior-posterior  (AP,  producing  a  stack  of  coronal  plane  images)  and  lateral -medial  (LM,  producing  a  stack  of  sagittal  plane  images)  views 
while  aligned  axially  (0°),  transversely  (90°)  or  obliquely  (23°)  to  the  superior-inferior  axis  of  the  vertebrae  (a  total  of  6  scans  per  spine).  The  vertebrae  were 
immersed  in  water  during  scans  to  simulate  soft  tissue. 

The  central  slices  from  AP-transverse  and  LM-transverse  (corresponding  to  coronal  and  sagittal  views)  were  analyzed  for  the  vertebral  heights  and  widths 
(Fig  1).  A  cuboidal  volume  of  interest  (VOI)  was  cropped  out  of  each  reconstructed  DTS  image-stacks  to  include  as  much  cancellous  bone  as  possible 
(Fig  .2).  Each  cropped  region  was  analyzed  using  the  fractal,  mean  intercept  length  (MIL)  and  line  fraction  deviation  (LFD)  methods  as  previously  described 
[2].  In  fractal  analysis,  fractal  dimension  (FD),  mean  lacunarity  (>i)  and  the  slope  of  lacunarity  vs.  box  size  relationship  (S^)  were  calculated  [5,  6].  In  MIL 
and  LFD  analyses,  the  mean  (Av),  standard  deviation  (SD),  maximum  (Max)  and  degree  of  anisotropy  (DA,  ratio  of  the  principal  measurements  of  MIL  or 
LFD)  were  recorded  for  each  slice  (Fig  3). 

Stepwise  forward  regression  method  was  used  to  construct  mixed  multiple  regression  models  that  examined  the  relationship  of  fractal,  MIL  and  LFD 
parameters  with  stiffness.  In  all  regression  models,  a  subject  variable  was  included  as  a  random  effect  to  account  for  pseudo-replication.  Separate  analyses 
were  performed  to  examine  contribution  of  DTS  derived  cancellous  bone  parameters  and  geometry  to  stiffness  with  and  without  the  presence  of  bone  mass 
information.  Each  scan  configuration  was  examined  separately.  Models  with  high  multicolinearity  (Variance  Inflation  Factor  >  5)  were  rejected. 

Results:  There  were  at  least  two  DTS  derived  cancellous  bone  parameters  for  each  scan  configuration  that  had  a  univariate  association  with  vertebral 
stiffness,  stronger  (R2adj  =  0.68-0.79)  than  that  between  stiffness  and  BV/TV  (R2adj  =  0.65).  Explained  variability  in  stiffness  was  improved  when  cancellous 
bone  parameters  were  used  in  multiple  regression  models  (Table  1),  with  the  axial  scans  in  the  LM  view  providing  highest  R2adj.  Inclusion  of  BV/TV  in  the 
cancellous-only  models  improved  R2adj  for  the  models  from  transverse  scans  but  not  from  other  orientations.  Addition  of  geometric  information  to 
cancellous-only  models  slightly  improved  the  models  from  the  oblique  and  axial  scans  but  not  from  the  transverse  scans.  The  most  explanatory  DTS -only 
model  (cancellous  bone  and  vertebral  geometry)  was  from  axial  scans  in  the  LM  view,  reaching  an  R2adj  value  of  0.90.  When  both  cancellous  and  geometry 
information  were  available  from  DTS,  addition  of  BV/TV  improved  the  models  for  transverse  scans  in  the  AP  view  only.  In  general,  a  measure  of  vertebral 
size  (LM  or  AP  width),  a  measure  of  trabecular  anisotropy  (MIL.DA  or,  similarly,  MIL.SD)  and  a  measure  of  trabecular  heterogeneity  (inter-planar  SD 
parameters),  all  measured  from  DTS  were  present  in  models  with  highest  R2adj. 

Discussion:  We  have  demonstrated  that  DTS  derived  cancellous  bone  microstructural  and  vertebral  geometry  parameters  are  strongly  associated  with  whole 
vertebral  stiffness.  DTS  alone  was  more  explanatory  for  vertebral  stiffness  than  BV/TV.  An  axial  scan  in  the  LM  view  appears  to  be  the  optimum 
configuration  for  a  model  of  vertebral  stiffness  with  the  highest  explanatory  capability,  if  DTS  is  to  be  used  alone.  Addition  of  BV/TV  did  not  improve  the 
explanatory  capability  of  the  models  in  this  case,  suggesting  that  additional  bone  mass  information  may  not  be  necessary  when  assessment  of  a  whole 
vertebra  is  made  using  DTS  in  the  LM-axial  configuration. 

Previous  studies  reported  an  association  of  cancellous  bone  stiffness  and  DTS  derived  cancellous  bone  microstructural  parameters  for  the  AP-transverse 
configuration  [1,2].  BV/TV  is  a  strong  predictor  of  cancellous  bone  stiffness  and  was  present  in  the  multiple  regression  models  for  cancellous  bone  stiffness. 
Consistent  with  the  previous  results  from  cancellous  bone,  models  of  whole  vertebral  stiffness  from  the  AP-transverse  configuration  reached  their  highest 
R2adj  when  BV/TV  was  included  in  the  model.  It  is  possible  then  that  cancellous  bone  can  be  assessed  using  DTS  alone  if  scan  configurations  other  than  the 
AP-transverse  are  considered. 

Previous  results  from  laboratory-scale  high  resolution  imaging  studies  (pCT)  suggested  that  cancellous  bone  microstructural  heterogeneity  is  associated  with 
cancellous  bone  stiffness,  independent  from  average  measures  of  micro  structure  [7].  By  using  the  entire  image  stack,  we  were  able  to  calculate  parameters 
that  represent  the  3D  heterogeneity  of  the  micro  structure.  Most  explanatory  models  of  stiffness  had  one  or  more  of  these  heterogeneity  parameters, 
demonstrating  a  major  advantage  of  DTS  over  imaging  methods  that  quantify  average  properties. 

The  predictor  variables  were  different  in  models  from  different  scan  configurations,  suggesting  that  complementary  information  may  be  obtained  from 
different  scan  configurations.  This  may  be  due  to  a  sensitivity  of  DTS  to  the  orientation  of  the  scanned  features  relative  to  the  orientation  of  the  scans  [8]. 
This  feature  may  be  developed  in  to  a  method  to  analyze  the  micro  structure  in  selected  orientations  with  DTS. 

Significance:  DTS  generated  images  analyzed  by  fractal,  MIL  and  LFD  methods  provided  good  to  excellent  predictions  of  in  vitro  vertebral  whole  bone 
stiffness.  This  can  be  further  improved  by  the  inclusion  of  vertebral  width  data  in  the  model. 
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Figure  1  -  Central  slice  from  anterior-posterior  (AP)  coronal  and  lateral-medial  (LM)  sagittal  scans  were  used  to  measure  waist-to-waist  width  and  endplate- 
to-endplate  height. 


Figure  2  -  Illustration  of  AP  and  LM  VOIs  in  a  vertebra.  The  red  and  the  blue  lines  represent  the  image  stacks  for  the  LM  and  the  AP  view,  respectively. 


Figure  3  -  [A]  A  square  region  cropped  out  of  the  DTS  image  is  analyzed  using  the  fractal,  MIL  and  LFD  analysis  methods.  MIL  and  LFD  are  calculated  in 
all  directions  in  the  image.  [B]  A  typical  polar  plot  representation  of  the  LFD  data,  from  which  an  average,  standard  deviation,  max  and  anisotropy  (LFD  at 


90°  /  LFD  at  0°)  can  be  calculated  for  each  slice.  [C]  Each  slice  in  the  stack  is  analyzed.  Stack  average  and  standard  deviation  (e.g.  SD  (LFD.SD))  are 
calculated  as  measures  of  3D  microstructure  within  the  bone. 


Table  1  -  Summary  of  the  most  explanatory  mixed  multiple  regression  models  for  vertebral  stiffness  for  each  scan  and  analysis  combination.  Each  cell 
shows  R2adj  followed  by  the  predictor(s)  in  the  model  (p-value).  Only  those  with  p  <  0.05  are  shown.  BV/TV  is  not  shown  for  clarity. 


DTS  cancellous  only 

DTS  only 

DTS  +  BV/TV 

DTS  +  BV/TV 

(cancellous  + 

(without  geometry) 

(with  geometry) 

geometry) 

AP 

Oblique 

0.81 

0.84 

0.79 

0.84 

(Anterior- 

posterior 

(23°) 

Av(FD)  (<0.02) 

AV(FD)  (0.04) 

Av(FD)  (0.004) 

AP  Width  (0.0001) 

coronal  plane) 

SD(LFD.SD)  (<0.0002) 

SD(LFD.SD)  (0.0001) 

SD(LFD.SD)  (0.0009) 

SD  LFD.SD  (0.006) 

Av(MIL.SD)  (<0.0004) 

AP  Height(<0.002) 

Av(MIL.SD)  (0.02) 

SD(MIL.Av)  (<0.02) 

Av(MIL.Av)  (0.002) 

Av(MIL.DA)  (<0.04) 

SD(MIL.Max)  (0.004) 

0.79 

0.79 

0.78 

0.82 

Axial 

(0°) 

SD(S?,)  (0.0001) 

SD(S?,)  (0.0001) 

SD(S^)  (0.0001) 

AP  Width  (0.007) 
SD(S0  (0.05) 

Transverse 

0.82 

0.82 

0.89 

0.91 

(90°) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.002) 

SD(MIL.DA)  (0.0002) 

SD(MIL.DA)  (0.0002) 

SD(MIL.DA)  (0.0001) 

SD(MIL.DA)  (0.0001) 

Av(MIL.Av)  (0.002) 

Av(MIL.Av)  (0.002) 

Av(MIL.Av)  (0.0001) 

Av(MIL.Av)  (0.0004) 

SD(MIL.SD)  (0.03) 

SD(MIL.SD)  (0.03) 

SD(MIL.SD)  (0.004) 

AP  Width  (0.009) 
SD(MIL.SD)  (0.02) 

LM 

Oblique 

0.76 

0.79 

0.77 

0.85 

(Lateral-medial 
sagittal  plane) 

(23°) 

SD(LFD.SD)  (0.002) 

LM  Width  (0.03) 

SD(LFD.SD)  (0.0004) 

LM  Width  (0.0001) 

SD(LFD.SD)  (0.003) 

SD(LFD.SD)  (0.0003) 

Axial 

0.88 

0.90 

0.83 

0.83 

(0°) 

Av(FD)  (0.0002) 

Av(FD)  (0.0001) 

Mean  (SO  (0.004) 

LM  Width  (0.006) 

SD(LFD.SD)  (0.0001) 

SD(LFD.SD)  (0.0001) 

SD(LFD.SD)  (0.0004) 

SD(LFD.Max)  (0.007) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.0001) 

Av  (SO  (0.05) 

Av(MIL.DA)  (0.002) 

Av(MIL.DA)  (0.005) 

SD(MIL.Av)  (0.004) 

SD(MIL.Max)  (0.008) 
LM  Width  (0.007) 

Transverse 

0.70 

0.70 

0.75 

0.76 

(90°) 

Av(MIL.SD)  (0.02) 

Av(MIL.SD)  (0.02) 

Av(MIL.SD)  (0.006) 

LM  Width  (0.0001) 

DTS  cancellous  only 

DTS  only 

DTS  +  BV/TV 

DTS  +  BV/TV 

(cancellous  + 

(without  geometry) 

(with  geometry) 

geometry) 

AP 

Oblique 

0.81 

0.84 

0.79 

0.84 

(Anterior- 

posterior 

(23°) 

Av(FD)  (0.02) 

AV(FD)  (0.04) 

Av(FD)  (0.004) 

AP  Width  (0.0001) 

coronal  plane) 

SD(LFD.SD)  (<0.0002) 

SD(LFD.SD)  (0.0001) 

SD(LFD.SD)  (0.0009) 

SD  LFD.SD  (0.006) 

Av(MIL.SD)  (0.0004) 

AP  Height(0.002) 

Av(MIL.SD)  (0.02) 

SD(MIL.Av)  (0.02) 

Av(MIL.Av)  (0.002) 

Av(MIL.DA)  (0.04) 

SD(MIL.Max)  (0.004) 

0.79 

0.79 

0.78 

0.82 

Axial 

(0°) 

SD(Sx)  (0.0001) 

SD(Sx)  (0.0001) 

SD(Sx)  (0.0001) 

AP  Width  (0.007) 
SD(S0  (0.05) 

Transverse 

0.82 

0.82 

0.89 

0.91 

(90°) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.002) 

SD(MIL.DA)  (0.0002) 

SD(MIL.DA)  (0.0002) 

SD(MIL.DA)  (0.0001) 

SD(MIL.DA)  (0.0001) 

Av(MIL.Av)  (0.002) 

Av(MIL.Av)  (0.002) 

Av(MIL.Av)  (0.0001) 

Av(MIL.Av)  (0.0004) 

SD(MIL.SD)  (0.03) 

SD(MIL.SD)  (0.03) 

SD(MIL.SD)  (0.004) 

AP  Width  (0.009) 
SD(MIL.SD)  (0.02) 

LM 

Oblique 

0.76 

0.79 

0.77 

0.85 

(Lateral-medial 
sagittal  plane) 

(23°) 

SD(LFD.SD)  (0.002) 

LM  Width  (0.03) 

SD(LFD.SD)  (0.0004) 

LM  Width  (0.0001) 

SD(LFD.SD)  (0.003) 

SD(LFD.SD)  (0.0003) 

Axial 

0.88 

0.90 

0.83 

0.83 

(0°) 

Av(FD)  (0.0002) 

Av(FD)  (0.0001) 

Mean  (SO  (0.004) 

LM  Width  (0.006) 

SD(LFD.SD)  (0.0001) 

SD(LFD.SD)  (0.0001) 

SD(LFD.SD)  (0.0004) 

SD(LFD.Max)  (0.007) 

Av(MIL.SD)  (0.0001) 

Av(MIL.SD)  (0.0001) 

Av  (SO  (0.05) 

Av(MIL.DA)  (0.002) 

Av(MIL.DA)  (0.005) 

SD(MIL.Av)  (0.004) 

SD(MIL.Max)  (0.008) 
LM  Width  (0.007) 

Transverse 

0.70 

0.70 

0.75 

0.76 

(90°) 

Av(MIL.SD)  (0.02) 

Av(MIL.SD)  (0.02) 

Av(MIL.SD)  (0.006) 

LM  Width  (0.0001) 

